calculation.cpp

Go to the documentation of this file.

/***************************************************************************
                        calculation.cpp  -  description
                             -------------------
    begin                : Thu April 28 2005
    copyright            : (C) 2005-2006 by Ben Swerts
    email                : bswerts@users.sourceforge.net
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/






// C++ header files
#include <cassert>
#include <iostream>

// Qt header files
#include <qapplication.h>
#include <qdir.h>
#include <qdom.h>
#include <qfile.h>
#include <qmessagebox.h>
#include <qprocess.h>
#include <qtimer.h>
#include <qurloperator.h>
#include <qwidget.h>

// Xbrabo header files
#include "atomset.h"
#include "calculation.h"
#include "crdfactory.h"
#include "domutils.h"
#include "globalbase.h"
#include "paths.h"


Calculation::Calculation(AtomSet* atomset, QWidget* parent, const char* name) : QObject(parent, name),
atoms(atomset),
parentWidget(parent),
braboInput(QStringList()),
stockInput(QStringList()),
relaxInput(QStringList()),
basissetList(QStringList()),
startVectorFile(QString::null),
atdensFile(QString::null),
affUpdateFreq(0),
calculationType(0),
crystalType(0),
updateBraboInput(true),
updateStockInput(false),
updateRelaxInput(false),
doBasissetCheck(true),
copyAtdens(false),
calcRunning(false),
calcPaused(false),
calcName(QString::null),
calcDir(QString::null),
calcProcess(0),
stopRequested(false),
backupFrequency(0),
calcModified(false),
calcError(NoError),
calcContinuable(false)
{
  assert(atomset != 0);
}

Calculation::~Calculation()
{

}

bool Calculation::start()
{
  if(calcRunning)
  {
    if(!calcPaused)
    {
      //QMessageBox::warning(parentWidget, tr("Start Calculation") , tr("Only 1 calculation can be run at a time"));
      return false;
    }
    else
    {
      calcPaused = false;
      doNextStep();
      return true;
    }
  }
  assert(!braboInput.isEmpty());
  if(calculationType == GlobalBase::GeometryOptimization)
    assert(!relaxInput.isEmpty());
  if(doBasissetCheck && !checkBasissets())
    return false;

  if(!checkExecutables())
    return false;

  if(!makeDirCurrent(calcDir, tr("Start Calculation")))
    return false;

  if(!copyStartVector())
    return false;

  bool calcContinue = false;
  QFile fileCRD(calcName + ".crd");
  if(calculationType == GlobalBase::GeometryOptimization && calcContinuable && fileCRD.exists())
  {
    // the optimization can possibly be continued from the current or previous cycle
    calcContinue = QMessageBox::question(parentWidget, tr("Geometry optimization"), 
      tr("The geometry optimization can possibly be continued."), //\nYou can choose to try to continue from the current cycle\nor the previous one."), 
      tr("Current cycle"), tr("Restart")) == 0;
  }
  if(!calcContinue && CrdFactory::writeToFile(atoms, calcName + ".crd", calcXF) != CrdFactory::OK)
  {
    QMessageBox::warning(parentWidget, tr("Start Calculation"), tr("Unable to write the initial coordinates"));
    return false;
  }

  if(calculationType == GlobalBase::GeometryOptimization)
  {
    QFile optcrd(calcName + ".c00");
    if(optcrd.exists())
    {
      if(!optcrd.remove())
      {
        QMessageBox::warning(parentWidget, tr("Start calculation"), tr("The existing .c00 file could not be removed."));
        return false;
      }
    }
  }
  
  QDir dir(calcDir);
  QStringList files = dir.entryList(calcName + "_*"); // all files named molecule_* 
  QFile file;
  for(QStringList::Iterator it = files.begin(); it != files.end(); it++)
  {
    file.setName(*it);
    file.remove();
  }

  // NOTE: after BRABO, STOCK and (MAFF,RELAX) can execute in parallel
  // NOTE: Always add all steps and decide whether or not to run them in doNextStep
  calculationSteps.clear();
  calculationSteps.push_back(BRABO);
  calculationSteps.push_back(STOCK);
  calculationSteps.push_back(UPDATE);
  if(calculationType == GlobalBase::GeometryOptimization)
  {
    calculationSteps.push_back(MAFF);
    calculationSteps.push_back(CNVRTAFF);
    calculationSteps.push_back(RELAX);
  }
  // add this step last
  calculationSteps.push_back(NEW_CYCLE);

  stopRequested = false;  
  calcRunning = true;
  currentCycle = 0;
  calcError = NoError;

  // start the calculation by invoking the dispatcher
  doNextStep();
  return true;
}

bool Calculation::pause()
{
  // return if the calculation is not running
  if(!calcRunning)
    return false;

  // start if calculation is already paused
  if(calcPaused)
    return start();

  // return if only 1 step is performed
  if(calculationSteps.size() == 2) 
    return false;

  calcPaused = true;
  setModified();
  return true;
}

bool Calculation::stop()
{
  if(!calcRunning)
    return false;
  
  if(calcPaused) // for loading a paused calculation from disk
  {
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return true;
  }

  int result = QMessageBox::information(parentWidget, tr("Stop calculation"), tr("Please specify how the calculation should be stopped."), 
               tr("After the current step"), tr("Immediately"), tr("Cancel"), 0, 2);
  switch(result)
  {
    case 0: // stop the calculation after the current step has ended
            stopRequested = true; 
            setModified();
            break;
    case 1: // stop the calculation immediately
            if(calcProcess != 0)
              calcProcess->kill();
            calcRunning = false;
            if(calcError == NoError)
              calcError = ManualStop;
            doNextStep();
            break;
    case 2: // cancel stopping
            return false;
  }
  return true;
}

void Calculation::setModified(const bool status)
{
  if(status && !calcModified)
    emit modified();
  calcModified = status;
}

void Calculation::setBraboInput(const QStringList stdInput, const QStringList basissets, const QString startVector, const bool preferStartVector, const unsigned int startVectorSize1, const unsigned int startVectorSize2)
{
  if(calcRunning && (calculationType == GlobalBase::SinglePointEnergy || calculationType == GlobalBase::EnergyAndForces))
    return;

  if(!stdInput.isEmpty())
  {
    braboInput = stdInput;
    updateBraboInput = true;
  }

  if(basissets != basissetList && !basissets.isEmpty())
  {
    basissetList = basissets;
    doBasissetCheck = true;
  }

  // - if the calculation is not running, check the validity of any existing startvector
  //   with the right name if preferStartVector = true, else use the provided startvector
  //   (or none if none is specified)
  // - if the calculation already started, only update the startvector if the basisset
  //   changed. 
  qDebug("setBraboInput: calcRunning = %d, prefer = %d, size1 = %d, size2 = %d",calcRunning, preferStartVector, startVectorSize1, startVectorSize2);
  qDebug("               startVector = " + startVector + ", startVectorFile = " + startVectorFile);
  if(!calcRunning)
  {
    if(preferStartVector)
    {
      QFileInfo fi(calcDir + QDir::separator() + calcName + ".sta");
      if(fi.size() == startVectorSize1 || fi.size() == startVectorSize2) // will fail if no filename like this exists
        startVectorFile = calcDir + QDir::separator() + calcName + ".sta";
      else
      {
        fi.setFile(startVector);
        qDebug("size of " + startVector + " = %d", fi.size());
        if(fi.size() == startVectorSize1 || fi.size() == startVectorSize2) // will always return false if startVector is empty
          startVectorFile = startVector; 
        else
          startVectorFile = QString::null; // do not use a starting vector at the next step (initial step here)
      }
    }
    else
    {
      QFileInfo fi(startVector);
      if(fi.size() == startVectorSize1 || fi.size() == startVectorSize2)
        startVectorFile = startVector;
      else
        startVectorFile = QString::null;
    }
  }
  else if(doBasissetCheck)
  {
    QFileInfo fi(startVectorFile);
    if(fi.size() != startVectorSize1 && fi.size() != startVectorSize2)
    {
      // the calculation is running, the basis set has changed and the
      // current starting vector cannot be used anymore.
      // WARNING: this does not work when the calculation is doing the first BRABO 
      //          call without a starting vector and before the first iteration of
      //          the SCF has been completed (thus before a starting vector has been 
      //          written). In that case startVectorFile is still empty.
      fi.setFile(startVector);
      if(fi.size() == startVectorSize1 || fi.size() == startVectorSize2)
        startVectorFile = startVector; // use the given start vector 
      else
        startVectorFile = QString::null; // regenerate the starting vector    
    }
  }
  qDebug("chosen startvector = " + startVectorFile);
}

void Calculation::setStockInput(const QStringList stin, const QString atdens)
{
  if(!stin.isEmpty())
  {
    stockInput = stin;
    updateStockInput = true;
  }

  if(!atdens.isEmpty())
  {
    atdensFile = atdens;
    copyAtdens = true;
  }
}

void Calculation::setRelaxInput(const QStringList aff, const QStringList maff, const unsigned int updateFreq, const unsigned int maxSteps, const std::vector<unsigned int>& steps, const std::vector<double>& factors)
{
  if(!aff.isEmpty())
  {
    relaxInput = aff;
    updateRelaxInput = true;
    maffInput = maff;
  }
  affUpdateFreq = updateFreq;
  if(!steps.empty() && steps.size() == factors.size())
  {
    scaleSteps.clear();
    scaleSteps.reserve(steps.size());
    scaleSteps.assign(steps.begin(), steps.end());
    scaleFactors.clear();
    scaleFactors.reserve(factors.size());
    scaleFactors.assign(factors.begin(), factors.end());
  }
  calcMaxCycle = maxSteps;
}

void Calculation::setCalculationType(const unsigned int type, const unsigned int buur)
{
  if(!calcRunning)
  {
    calculationType = type;
    crystalType = buur;
  }
}

void Calculation::setParameters(const QString name, const QString dir, const bool format)
{
  calcName = name;
  calcDir = dir;
  calcXF = format;
}

void Calculation::setBackup(const unsigned int freq, const bool brabo, const bool stock, const bool relax, const bool aff, const bool crd)
{
  backupFrequency = freq;
  if(backupFrequency != 0)
  {
    backupBrabo = brabo;
    backupStock = stock;
    backupRelax = relax;
    backupAFF = aff;
    backupCRD = crd;
  }
  setModified();
}

void Calculation::setContinuable(const bool status)
{
  calcContinuable = status;
}

bool Calculation::isRunning() const
{
  return calcRunning;
}

bool Calculation::isPaused() const
{
  return calcPaused;
}

void Calculation::getRefinementParameters(bool& largestCart, bool& magnCart, bool& largestInt, bool& magnInt, bool& largestShift)
{
  largestCart = false;
  magnCart = false;
  largestInt = false;
  magnInt = false;
  largestShift = false;

  QFile file(calcDir + QDir::separator() + calcName + ".aou");
  if(!file.open(IO_ReadOnly))
    return;
  QTextStream stream(&file);
  while(!stream.atEnd())
  {
    QString line = stream.readLine();
    if(line.contains("OK"))
    {
      if(line.contains("LARGEST CARTESIAN FORCE"))
        largestCart = true;
      else if(line.contains("MAGNIT. CART. FORCE VEC"))
        magnCart = true;
      else if(line.contains("LARGEST INTERNAL FORCE"))
        largestInt = true;
      else if(line.contains("MAGNIT. INT. FORCE VEC"))
        magnInt = true;
      else if(line.contains("LARGEST SHIFT"))
        largestShift = true;
    }
  }
}

void Calculation::getAvailableOutputs(std::vector<unsigned int>& cycles, std::vector<bool>& brabo, std::vector<bool>& stock, std::vector<bool>& relax, std::vector<bool>& aff)
{
  cycles.clear();
  brabo.clear();
  stock.clear();
  relax.clear();
  aff.clear();

  bool out, stou, aou, aff2;
  QFile file;
  QString prefix = calcDir + QDir::separator() + calcName + "_";
  for(unsigned int i = 1; i <= currentCycle; i++)
  {
    QString prefixCycle = prefix + QString::number(i);
    file.setName(prefixCycle + ".out");
    out = file.exists();
    file.setName(prefixCycle + ".stou");
    stou = file.exists();
    file.setName(prefixCycle + ".aou");
    aou = file.exists();
    file.setName(prefixCycle + ".aff");
    aff2 = file.exists();
    if(out || stou || aou || aff2)
    {
      // at least one output exists so add a cycle
      cycles.push_back(i);
      brabo.push_back(out);
      stock.push_back(stou);
      relax.push_back(aou);
      aff.push_back(aff2);
    }
  }
}

QStringList Calculation::braboOutput(const unsigned int step)
{
  return output(".out", step);
}

QStringList Calculation::stockOutput(const unsigned int step)
{
  return output(".stou", step);
}

QStringList Calculation::relaxOutput(const unsigned int step)
{
  return output(".aou", step);
}

QStringList Calculation::affOutput(const unsigned int step)
{
  return output(".aff", step);
}

void Calculation::loadCML(const QDomElement* root)
{
  const QString prefix = "calculation_";
  QDomNode childNode = root->firstChild();
  while(!childNode.isNull())
  {
    if(childNode.isElement() && childNode.toElement().tagName() == "parameter")
    {
      if(DomUtils::dictEntry(childNode, prefix + "starting_vector"))
        DomUtils::readNode(&childNode, &startVectorFile);
      else if(DomUtils::dictEntry(childNode, prefix + "update_energy_and_forces"))
        DomUtils::readNode(&childNode, &updateBraboInput);
      else if(DomUtils::dictEntry(childNode, prefix + "update_stockholder"))
        DomUtils::readNode(&childNode, &updateStockInput);
      else if(DomUtils::dictEntry(childNode, prefix + "update_geometry_optimization"))
        DomUtils::readNode(&childNode, &updateRelaxInput);
      else if(DomUtils::dictEntry(childNode, prefix + "check_basissets"))
        DomUtils::readNode(&childNode, &doBasissetCheck);
      else if(DomUtils::dictEntry(childNode, prefix + "running"))
        DomUtils::readNode(&childNode, &calcRunning);
      else if(DomUtils::dictEntry(childNode, prefix + "paused"))
        DomUtils::readNode(&childNode, &calcPaused);
      else if(DomUtils::dictEntry(childNode, prefix + "current_cycle"))
        DomUtils::readNode(&childNode, &currentCycle);
      else if(DomUtils::dictEntry(childNode, prefix + "error"))
        DomUtils::readNode(&childNode, &calcError);
      else if(DomUtils::dictEntry(childNode, prefix + "continuable"))
        DomUtils::readNode(&childNode, &calcContinuable);
      else if(DomUtils::dictEntry(childNode, prefix + "steps"))
      {
        std::vector<unsigned int> steps;
        DomUtils::readNode(&childNode, &steps);
        calculationSteps.clear();
        for(unsigned int i = 0; i < steps.size(); i++)
        {
          switch(steps[i])
          {
            case NEW_CYCLE: calculationSteps.push_back(NEW_CYCLE);
                            break;
            case ACHAR:     calculationSteps.push_back(ACHAR);
                            break;
            case BRABO:     calculationSteps.push_back(BRABO);
                            break;
            case BUUR:      calculationSteps.push_back(BUUR);
                            break;
            case CNVRTAFF:  calculationSteps.push_back(CNVRTAFF);
                            break;
            case DISTOR:    calculationSteps.push_back(DISTOR);
                            break;
            case FORKON:    calculationSteps.push_back(FORKON);
                            break;
            case MAFF:      calculationSteps.push_back(MAFF);
                            break;
            case RELAX:     calculationSteps.push_back(RELAX);
                            break;
            case STOCK:     calculationSteps.push_back(STOCK);
                            break;
            case UPDATE:    calculationSteps.push_back(UPDATE);
                            break;
          }
        }
      }
    }
    childNode = childNode.nextSibling();
  }
}

void Calculation::saveCML(QDomElement* root)
{
  const QString prefix = "calculation_";
  DomUtils::makeNode(root, startVectorFile, prefix + "starting_vector");
  DomUtils::makeNode(root, updateBraboInput, prefix + "update_energy_and_forces");
  DomUtils::makeNode(root, updateStockInput, prefix + "update_stockholder");
  DomUtils::makeNode(root, updateRelaxInput, prefix + "update_geometry_optimization");
  DomUtils::makeNode(root, doBasissetCheck, prefix + "check_basissets");
  DomUtils::makeNode(root, calcRunning, prefix + "running");
  DomUtils::makeNode(root, calcPaused, prefix + "paused");
  DomUtils::makeNode(root, currentCycle, prefix + "current_cycle");
  DomUtils::makeNode(root, calcError, prefix + "error");
  DomUtils::makeNode(root, calcContinuable, prefix + "continuable");
  std::vector<unsigned int> calcSteps;
  for(std::list<CalculationStep>::iterator it = calculationSteps.begin(); it != calculationSteps.end(); it++)
    calcSteps.push_back(static_cast<unsigned int>(*it));
  DomUtils::makeNode(root, calcSteps, prefix + "steps");
}

void Calculation::writeInput()
{
  if(!makeDirCurrent(calcDir, tr("Write all input files")))
    return;

  QFile fileInp(calcName + ".inp");
  if(!fileInp.open(IO_WriteOnly))
  {
    QMessageBox::warning(parentWidget, tr("Write all input files"), tr("Unable to write the file") +
      "\n" + fileInp.name());
    return;
  }
  QTextStream streamInp(&fileInp);
  streamInp << QString(braboInput.join("\n")+"\n");
  fileInp.close();
  if(calcRunning)
    return; // when a calculation is running, all other files are already present

  if(CrdFactory::writeToFile(atoms, calcName + ".crd", calcXF) != CrdFactory::OK)
  {
    QMessageBox::warning(parentWidget, tr("Write all input files"), tr("Unable to write the file") + "\n" + calcName + ".crd");
    return;
  }

  if(!copyStartVector())
    return;

  if(!stockInput.isEmpty())
  {
    QFile file(calcName + ".stin");
    if(!file.open(IO_WriteOnly))
    {
      QMessageBox::warning(parentWidget, tr("Write all input files"), tr("Unable to write the file") +
      "\n" + file.name());
      return;
    }
    QTextStream stream(&file);
    stream << QString(stockInput.join("\n")+"\n");
    fileInp.close();
  }

  if(!maffInput.isEmpty())
  {
    QFile file(calcName + ".maffinp");
    if(!file.open(IO_WriteOnly))
    {
      QMessageBox::warning(parentWidget, tr("Write all input files"), tr("Unable to write the file") +
      "\n" + file.name());
      return;
    }
    QTextStream stream(&file);
    stream << QString(maffInput.join("\n")+"\n");
    fileInp.close();
  }

  if(!relaxInput.isEmpty())
  {
    QString extension = ".aff";
    if(!maffInput.isEmpty())
      extension += "1";
    QFile file(calcName + extension);
    if(!file.open(IO_WriteOnly))
    {
      QMessageBox::warning(parentWidget, tr("Write all input files"), tr("Unable to write the file") +
      "\n" + file.name());
      return;
    }
    QTextStream stream(&file);
    stream << QString(relaxInput.join("\n")+"\n");
    fileInp.close();
  }
}

void Calculation::clean()
{
  if(calcRunning)
    return;

  QDir dir(calcDir);
  if(!dir.exists())
    return; // certainly clean

  if(!makeDirCurrent(calcDir, tr("Clean calculation")))
    return;

  QStringList files = dir.entryList(calcName + ".*"); // all in and output files
  files += dir.entryList(calcName + "_*"); // all saved outputs
  files += dir.entryList("fort.*"); // all fortran unnnamed unit files
  QFile file;
  for(QStringList::Iterator it = files.begin(); it != files.end(); it++)
  {
    file.setName(*it);
    file.remove();
  }
  dir.setPath(calcDir + QDir::separator() + "tmp");
  if(!dir.exists())
    return;
  dir.setCurrent(dir.path());
  files = dir.entryList(calcName + ".*"); // all in and output files
  files += dir.entryList("fort.*"); // all fortran unnnamed unit files
  for(QStringList::Iterator it = files.begin(); it != files.end(); it++)
  {
    file.setName(*it);
    file.remove();
  }  
  dir.cdUp();
  dir.rmdir(calcDir + QDir::separator() + "tmp"); // will only be removed if it is empty 
}


void Calculation::finishBrabo()
{
  assert(calcProcess != 0);

  if(!calcRunning)
    return;

  if(calcProcess->isRunning())
  {
    calcProcess->tryTerminate();
    QTimer::singleShot(1000, this, SLOT(finishBrabo));
    return;
  }
  
  bool noError = calcProcess->normalExit();
  delete calcProcess;
  calcProcess = 0;

  if(!noError)
  {
    if(calculationType == GlobalBase::SinglePointEnergy)
      QMessageBox::warning(parentWidget, tr("Energy"), tr("A problem occured during the calculation of the energy"));
    else
      QMessageBox::warning(parentWidget, tr("Energy & Forces"), tr("A problem occured during the calculation of the energy & forces"));    
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
  }
  
  doNextStep();
}

void Calculation::readStdOutBrabo()
{
  assert(calcProcess != 0);

  if(!calcRunning)
    return;
  
  while(calcProcess->canReadLineStdout())
  {
    QString line = calcProcess->readLineStdout();
    //qDebug("stdout: " + line);
    if(line.length() == 75 && line.left(30) != "                              ")
    {
      bool oki, oke;
      unsigned int iteration = line.left(4).toUInt(&oki);
      double energy = line.mid(4,16).toDouble(&oke);
      if(oki && oke)
        emit newIteration(iteration, energy);
      if(oki && oke && startVectorFile.isEmpty())
        startVectorFile = calcDir + QDir::separator() + calcName + ".sta"; // makes sure this variable is filled as soon as the actual file exists
    }
    if(calcError == NoError && line.left(58) == "                                        NO CONVERGENCE !!!")
      calcError = NoConvergence; // the only error that can be deducted from stdout
  }
  
/*      // keep the first error that is encountered
      if(line == "        ERROR DETECTED        ERROR DETECTED        ERROR DETECTED        ERROR DETECTED        ERROR DETECTED")
        calcError = UndefinedError;
      else if(line.left(7) == " NUCLEI" && line.mid(17,10) == "TOO CLOSE,")
        calcError = CloseNuclei;
      qDebug("left7 = |" + line.left(7) + "|, mid = |" + line.mid(17,10) + "|");
*/
}

void Calculation::finishStock()
{
  assert(calcProcess != 0);

  if(!calcRunning)
    return;

  if(calcProcess->isRunning())
  {
    calcProcess->tryTerminate();
    QTimer::singleShot(1000, this, SLOT(finishStock));
    return;
  }
  
  bool noError = calcProcess->normalExit();
  delete calcProcess;
  calcProcess = 0;

  if(!noError)
  {
    QMessageBox::warning(parentWidget, tr("Calculation of stockholder charges"), tr("A problem occured during the calculation of the stockholder charges"));    
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }

  QFile inputFile(calcDir + QDir::separator() + "fort.7");
  QFile appendFile(calcDir + QDir::separator() + calcName + ".pun");
  if(!inputFile.open(IO_ReadOnly) || !appendFile.open(IO_WriteOnly | IO_Append))
  {
    QMessageBox::warning(parentWidget, tr("Calculation of stockholder charges"), tr("A problem occured during the calculation of the stockholder charges"));    
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }
  appendFile.writeBlock(inputFile.readAll());
  inputFile.close();
  appendFile.close();

  doNextStep();
}

void Calculation::finishMAFF()
{
  assert(calcProcess != 0);

  if(!calcRunning)
    return;

  if(calcProcess->isRunning())
  {
    calcProcess->tryTerminate();
    QTimer::singleShot(1000, this, SLOT(finishMAFF));
    return;
  }
  
  bool noError = calcProcess->normalExit();
  delete calcProcess;
  calcProcess = 0;

  if(!noError)
  {
    QMessageBox::warning(parentWidget, tr("Updating the internal coordinates"), tr("A problem occured during the generation of new internal coordinates"));    
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
  }

  doNextStep();
}

void Calculation::finishCnvrtAFF()
{
  assert(calcProcess != 0);

  if(!calcRunning)
    return;

  if(calcProcess->isRunning())
  {
    calcProcess->tryTerminate();
    QTimer::singleShot(1000, this, SLOT(finishCnvrtAFF));
    return;
  }
  
  bool noError = calcProcess->normalExit();
  delete calcProcess;
  calcProcess = 0;

  if(!noError)
  {
    QMessageBox::warning(parentWidget, tr("Updating the internal coordinates"), tr("A problem occured during the conversion of the new internal coordinates"));    
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }

  QFile fileNewAFF(calcDir + QDir::separator() + calcName + ".aff_new");
  QFile fileAFF(calcDir + QDir::separator() + calcName + ".aff");
  if(!fileNewAFF.open(IO_ReadOnly) || !fileAFF.open(IO_WriteOnly | IO_Translate)) // suddenly Relax wants files with crlf on Windows?
  {
    QMessageBox::warning(parentWidget, tr("Updating the internal coordinates"), tr("A problem occured during combining the new internal coordinates"));    
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }
  QTextStream streamNewAFF(&fileNewAFF);
  QTextStream streamAFF(&fileAFF);
  streamAFF << QString(relaxInput.join("\n")+"\n");
  QString line;
  while(!line.contains("GBMA", false) && !line.contains("BMAT", false))
    line = streamNewAFF.readLine();
  streamAFF << streamNewAFF.read(); // add the remainder of the file
  fileNewAFF.close();
  fileAFF.close();

  doNextStep();
}

void Calculation::finishRelax()
{
  assert(calcProcess != 0);

  if(!calcRunning)
    return;

  if(calcProcess->isRunning())
  {
    calcProcess->tryTerminate();
    QTimer::singleShot(1000, this, SLOT(finishRelax));
    return;
  }
  
  bool noError = calcProcess->normalExit();
  delete calcProcess;
  calcProcess = 0;

  qDebug("finishRelax: noError = %d", noError);
  if(!noError)
  {
    QMessageBox::warning(parentWidget, tr("Geometry optimization"), tr("A problem occured during the calculation of new coordinates"));    
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }

  QFile fileC00(calcDir + QDir::separator() + calcName + ".c00");
  if(fileC00.exists())
  {
    // read the coordinates
    unsigned int result = CrdFactory::readFromFile(atoms, fileC00.name());
    if(result != CrdFactory::OK)
    {
      QMessageBox::warning(parentWidget, tr("Run calculation"), tr("A problem occured reading the optimized coordinates."));
      if(calcError == NoError)
        calcError = UndefinedError;
    }
    else
      calcError = NoError;
    calcRunning = false;
    doNextStep();
    return;
  }

  QFile fileNCR(calcDir + QDir::separator() + calcName + ".ncr");
  if(!fileNCR.exists())
  {
    QMessageBox::warning(parentWidget, tr("Run calculation"), tr("A problem occured reading the new coordinates\n(the coordinate file does not exist)."));
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }

  QFile fileCrd(calcDir + QDir::separator() + calcName + ".crd");
  if(!fileCrd.remove()) // this is suddenly needed otherwise the rename won't work on Windows
  {
    QMessageBox::warning(parentWidget, tr("Run calculation"), tr("A problem occured updating the coordinates."));
    calcRunning  = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }
  QUrlOperator local(calcDir);
  const QNetworkOperation* operation = local.rename(calcName + ".ncr", calcName + ".crd");
  while((operation->state() == QNetworkProtocol::StWaiting) || (operation->state() == QNetworkProtocol::StInProgress))
    qApp->processEvents(); // shouldn't do much as renaming should be very fast as opposed to copying
  if(operation->state() != QNetworkProtocol::StDone)
  {
    qDebug("operation->state() = %d",operation->state());
    qDebug(" error code = %d",operation->errorCode());
    qDebug("error string = "+operation->protocolDetail());
    QMessageBox::warning(parentWidget, tr("Run calculation"), tr("A problem occured updating the coordinates."));
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
  }

  doNextStep();
}

void Calculation::readStdOutAll()
{
  assert(calcProcess != 0);

  if(!calcRunning)
    return;
  
  while(calcProcess->canReadLineStdout())
    calcProcess->readLineStdout();
}


bool Calculation::checkBasissets()
{
  QStringList missingFiles;
  for(QStringList::iterator it = basissetList.begin(); it != basissetList.end(); it++)
  {
    QFile file(*it);
    if(!file.exists())
      missingFiles += *it;
  }
  if(missingFiles.isEmpty())
  {
    doBasissetCheck = false;
    return true;
  }
  
  QMessageBox::warning(parentWidget, tr("Check basis sets") , tr("The following basis set files could not be found:") 
                       + "\n\n - " + missingFiles.join("\n - "));
  return false;
}

bool Calculation::checkExecutables()
{
  QStringList missingExecutables;
  QFile file(Paths::brabo);
  if(!file.exists())
    missingExecutables += "Brabo (" + Paths::brabo  +")";
  if(!stockInput.isEmpty())
  {
    file.setName(Paths::relax);
    if(!file.exists())
      missingExecutables += "Stock (" + Paths::stock  +")";
  }
  if(calculationType == GlobalBase::GeometryOptimization)
  {
    file.setName(Paths::relax);
    if(!file.exists())
      missingExecutables += "Relax (" + Paths::relax  +")";
    if(!maffInput.isEmpty())
    {
      file.setName(Paths::maff);
      if(!file.exists())
        missingExecutables += "Maff (" + Paths::maff  +")";
      file.setName(Paths::cnvrtaff);
      if(!file.exists())
        missingExecutables += "Cnvrtaff (" + Paths::cnvrtaff +")";
    }
  }
  if(missingExecutables.isEmpty())
    return true;

  QMessageBox::warning(parentWidget, tr("Check executables") , tr("The following executable files could not be found:") 
                       + "\n\n - " + missingExecutables.join("\n - ") + "\n\n" + tr("Please check the Preferences."));
  return false;
}

bool Calculation::makeDirCurrent(const QString dir, const QString title)
{
  QDir workDir = dir;
  if(!workDir.exists())
  {
    if(!workDir.mkdir(workDir.path()))
    {
      QMessageBox::warning(parentWidget, title, tr("Unable to create the directory ") + workDir.path());
      return false;
    }
  }
  if(!QDir::setCurrent(workDir.path()))
  {
    QMessageBox::warning(parentWidget, title, tr("Unable to change to the directory ") + workDir.path());
    return false;
  }
  return true;
}

bool Calculation::copyStartVector()
{
  if(startVectorFile.isEmpty())
    return true;

  if(startVectorFile == calcDir + QDir::separator() + calcName + ".sta")
    return true;

  qDebug("copying startvector");
  if(!copyFile(startVectorFile, calcDir + QDir::separator() + calcName + ".sta"))
  {
    QMessageBox::warning(parentWidget, tr("Copy starting vector"), tr("Unable to copy the starting vector \n") 
                         + startVectorFile + tr(" to the calculation directory"));
    return false;
  }
  return true;
}

void Calculation::doNextStep()
{
  qDebug("doNextStep: calcPaused = %d, calcRunning = %d, stopRequested = %d",
                      calcPaused, calcRunning, stopRequested);

  // do nothing if the calculation is paused
  if(calcPaused)
    return;

  setModified();

  // check whether the calculation has finished successfully or an error occured
  if(!calcRunning || stopRequested)
  {
    runCleanup();
    return;
  }

  // get the current step
  CalculationStep currentStep = calculationSteps.front();
  // take it from the front and insert it at the back for the next call of doNextStep
  calculationSteps.pop_front();
  calculationSteps.push_back(currentStep);

  // do the appropriate action 
  const unsigned int thisStep = currentStep;
  switch(thisStep) // switch doesn't work with enums
  {
    case NEW_CYCLE: // start of a new cycle
                    qDebug("doNextStep: currentStep = NEW_CYCLE");
                    if(calculationType == GlobalBase::SinglePointEnergy || calculationType == GlobalBase::EnergyAndForces)
                      calcRunning = false;
                    else if(currentCycle == calcMaxCycle)
                    {
                      calcRunning = false;
                      if(calcError == NoError)
                        calcError = MaxCyclesExceeded;
                    }
                    else
                    {
                      // backup the output files
                      backupOutputs();
                      // start a new optimization cycle
                      currentCycle++;
                      emit newCycle(currentCycle+1);
                      doNextStep();
                      return;
                    }
                    break;

    case BRABO:     // run BRABO     
                    qDebug("doNextStep: currentStep = BRABO");
                    runBrabo();
                    break;

    case CNVRTAFF:  // run CNVRTAFF
                    qDebug("doNextStep: currentStep = CNVRTAFF");
                    if(!maffInput.isEmpty() && ((affUpdateFreq == 0 && currentCycle == 0) || (affUpdateFreq != 0 && (currentCycle % affUpdateFreq == 0) || updateRelaxInput)))
                      // only run CnvrtAFF/MAFF if
                      // - there is an input file for maff (if not, no automatic generation of IC's is requested
                      // - it's the first cyle and only 1 AFF file is to be made (when affUpdateFreq == 0)
                      // - the AFF file should be regenerated every affUpdateFreq cycles
                      // - automatic regeneration and a change in the AFF header file
                      runCnvrtAFF();
                    else
                    {
                      doNextStep();
                      return;
                    }
                    break;

    case MAFF:      // run MAFF
                    qDebug("doNextStep: currentStep = MAFF");
                    if(!maffInput.isEmpty() && ((affUpdateFreq == 0 && currentCycle == 0) || (affUpdateFreq != 0 && (currentCycle % affUpdateFreq == 0) || updateRelaxInput)))
                      // only run CnvrtAFF/MAFF if
                      // - there is an input file for maff (if not, no automatic generation of IC's is requested
                      // - it's the first cyle and only 1 AFF file is to be made (when affUpdateFreq == 0)
                      // - the AFF file should be regenerated every affUpdateFreq cycles
                      // - automatic generation and a change in the AFF header file
                      runMAFF();
                    else
                    {
                      doNextStep();
                      return;
                    }
                    break;

    case RELAX:     // run RELAX
                    qDebug("doNextStep: currentStep = RELAX");
                    runRelax();
                    break;
         
    case STOCK:     // run STOCK
                    qDebug("doNextStep: currentStep = STOCK");
                    if(stockInput.isEmpty())
                    {
                      doNextStep();
                      return;
                    }
                    else
                      runStock(); // only run it when these charges are requested
                    break;

    case UPDATE:    // run UPDATE
                    qDebug("doNextStep: currentStep = UPDATE");
                    runUpdate();
                    return; // never call runCleanup after runUpdate (runCleanup never changes calcRunning)
                    break;
  }
  // do a cleanup if something went wrong or the end of the calculation, 
  // signalled by setting calcRunning to false
  if(!calcRunning)
    runCleanup();
}

void Calculation::runBrabo()
{
  assert(calcProcess == 0);

  if(updateBraboInput)
  {
    // update this input if it should exist as a physical file on disk

    updateBraboInput = false;
  }

  QString stdInput;
  if(startVectorFile.isEmpty())
  {
    for(QStringList::Iterator it = braboInput.begin(); it != braboInput.end(); it++)
    {
      if((*it).left(4) != "star")
        stdInput += *it + "\n";
    }
  }
  else
    stdInput = braboInput.join("\n")+"\n";

  if(!copyStartVector())
  {
    calcRunning = false;
    if(calcError == NoError)
      calcError = UndefinedError;
    doNextStep();
    return;
  }

  calcProcess = new QProcess(this);
  calcProcess->setWorkingDirectory(calcDir);
  calcProcess->addArgument(Paths::brabo);
  connect(calcProcess, SIGNAL(processExited()), this, SLOT(finishBrabo()));
  connect(calcProcess, SIGNAL(readyReadStdout()), this, SLOT(readStdOutBrabo()));
  calcProcess->launch(stdInput);
}

void Calculation::runStock()
{
  assert(calcProcess == 0);

  if(updateStockInput)
  {
    QFile file(calcDir + QDir::separator() + calcName + ".stin");
    if(!file.open(IO_WriteOnly))
    {
      QMessageBox::warning(parentWidget, tr("Start the calculation of stockholder charges"), tr("Unable to update the file") +
        "\n" + file.name());
      calcRunning = false;
      if(calcError == NoError)
        calcError = UndefinedError;
      doNextStep();
      return;
    }
    QTextStream stream(&file);
    stream << QString(stockInput.join("\n")+"\n");
    file.close();

    if(!atdensFile.isEmpty())
    {
      file.setName(calcDir + QDir::separator() + calcName + ".atdens");
      if(!file.open(IO_WriteOnly | IO_Translate))
      {
        QMessageBox::warning(parentWidget, tr("Start the calculation of stockholder charges"), tr("Unable to update the file") +
          "\n" + file.name());
        calcRunning = false;
        if(calcError == NoError)
          calcError = UndefinedError;
        doNextStep();
        return;
      }
      stream.setDevice(&file);
      stream << atdensFile;
      file.close();
      atdensFile = QString::null;
    }
    updateStockInput = false;
  }

  calcProcess = new QProcess(this);
  calcProcess->setWorkingDirectory(calcDir);
  calcProcess->addArgument(Paths::stock);
  connect(calcProcess, SIGNAL(processExited()), this, SLOT(finishStock()));
  connect(calcProcess, SIGNAL(readyReadStdout()), this, SLOT(readStdOutAll()));
  calcProcess->launch(calcName + "\n");
}

void Calculation::runMAFF()
{
  assert(calcProcess == 0);
  assert(!maffInput.isEmpty());

  calcProcess = new QProcess(this);
  calcProcess->setWorkingDirectory(calcDir);
  calcProcess->addArgument(Paths::maff);
  connect(calcProcess, SIGNAL(processExited()), this, SLOT(finishMAFF()));
  connect(calcProcess, SIGNAL(readyReadStdout()), this, SLOT(readStdOutAll()));
  calcProcess->launch(maffInput.join("\n")+"\n"); 
}

void Calculation::runCnvrtAFF()
{
  assert(calcProcess == 0);

  updateRelaxInput = false; // should be disabled here because if in runMAFF this 
                            // routine wouldn't be called and if not at all
                            // an update in runRelax would overwrite stuff.

  calcProcess = new QProcess(this);
  calcProcess->setWorkingDirectory(calcDir);
  calcProcess->addArgument(Paths::cnvrtaff);
  connect(calcProcess, SIGNAL(processExited()), this, SLOT(finishCnvrtAFF()));
  connect(calcProcess, SIGNAL(readyReadStdout()), this, SLOT(readStdOutAll()));
  calcProcess->launch(calcName + "\n");
}

void Calculation::runRelax()
{
  assert(calcProcess == 0);

  if(updateRelaxInput)
  {
    QFile file(calcDir + QDir::separator() + calcName + ".aff");
    if(!file.open(IO_WriteOnly))
    {
      QMessageBox::warning(parentWidget, tr("Updating the internal coordinates"), tr("A problem occured writing the new internal coordinates"));    
      calcRunning = false;
      if(calcError == NoError)
        calcError = UndefinedError;
      doNextStep();
      return;
    }
    QTextStream stream(&file);
    stream << QString(relaxInput.join("\n") + "\n");
    file.close();    
    updateRelaxInput = false;
  }

  calcProcess = new QProcess(this);
  calcProcess->setWorkingDirectory(calcDir);
  calcProcess->addArgument(Paths::relax);
  connect(calcProcess, SIGNAL(processExited()), this, SLOT(finishRelax()));
  connect(calcProcess, SIGNAL(readyReadStdout()), this, SLOT(readStdOutAll()));
  calcProcess->launch(calcName + "\n" + QString::number(currentCycle+1) + "\nn\n" + QString::number(scaleFactor(currentCycle)) +"\n");
}

void Calculation::runCleanup()
{
  qDebug("Calculation::runCleanup() called");

  if(calcError == NoError)
    setContinuable(false); // when the structure is optimized, it should never be continued

  emit finished(calcError);
}

void Calculation::runUpdate()
{
  if(calculationType == GlobalBase::GeometryOptimization && currentCycle != 0)
  {
    QFile file(calcDir + QDir::separator() + calcName + ".crd");
    if(!file.exists() || CrdFactory::readFromFile(atoms, file.name()) != CrdFactory::OK)
      QMessageBox::warning(parentWidget, tr("Run calculation"), tr("The view could not be updated with the new coordinates."));
  }

  if(calculationType != GlobalBase::SinglePointEnergy)
  {
    QFile file(calcDir + QDir::separator() + calcName + ".pun");    
    if(!file.exists() || CrdFactory::readForces(atoms, file.name()) != CrdFactory::OK)
    {
      QMessageBox::warning(parentWidget, tr("Run calculation"), tr("The view could not be updated with the new forces."));
      qDebug("filename: "+file.name());
      qDebug("return value: %d", CrdFactory::readForces(atoms, file.name())); 
    }
  }

  atoms->setCharges(loadFromPunch("MULL", atoms->count(), 10, 8), AtomSet::Mulliken);

  atoms->setCharges(loadFromPunch("STOC", atoms->count(), 10, 8), AtomSet::Stockholder);

  emit updated();
  doNextStep();
}

double Calculation::scaleFactor(const unsigned int step)
{
  double result;
  unsigned int iterStep = 0;
  result = scaleFactors[0];
  for(unsigned int i = 0; i < scaleSteps.size(); i++)
  {
    iterStep += scaleSteps[i];
    if(step <= iterStep)
      result = scaleFactors[i];
    else
      break;
  }
  qDebug("Calculation::scaleFactor returns %f", result);
  for(unsigned int i = 0; i < scaleSteps.size(); i++)
    qDebug(" %d: step %d with factor %f",i,scaleSteps[i],scaleFactors[i]);

  return result;
}

void Calculation::backupOutputs()
{
  if(backupFrequency == 0 || currentCycle % backupFrequency != 0)
    return;

  if(backupBrabo)
    copyFile(calcDir + QDir::separator() + calcName + ".out", 
             calcDir + QDir::separator() + calcName + "_" + QString::number(currentCycle+1) + ".out");
  if(backupStock && !stockInput.isEmpty())
    copyFile(calcDir + QDir::separator() + calcName + ".stou", 
             calcDir + QDir::separator() + calcName + "_" + QString::number(currentCycle+1) + ".stou");
  if(backupRelax && calculationType == GlobalBase::GeometryOptimization)
    copyFile(calcDir + QDir::separator() + calcName + ".aou", 
             calcDir + QDir::separator() + calcName + "_" + QString::number(currentCycle+1) + ".aou");
  if(backupAFF && calculationType == GlobalBase::GeometryOptimization)
    copyFile(calcDir + QDir::separator() + calcName + ".aff", 
             calcDir + QDir::separator() + calcName + "_" + QString::number(currentCycle+1) + ".aff");
  if(backupCRD)
    copyFile(calcDir + QDir::separator() + calcName + ".crd", 
             calcDir + QDir::separator() + calcName + "_" + QString::number(currentCycle+1) + ".crd");
}

bool Calculation::copyFile(const QString source, const QString destination)
{
  QFile inputFile(source);
  QFile outputFile(destination);
  if(!inputFile.open(IO_ReadOnly) || !outputFile.open(IO_WriteOnly))
    return false;
  
  outputFile.writeBlock(inputFile.readAll());
  inputFile.close();
  outputFile.close();
  return true;
}

QStringList Calculation::output(const QString extension, const unsigned int step)
{
  QStringList result;
  QString fileName;
  if(step == 0)
    fileName = calcName + extension;
  else
    fileName = calcName + "_" + QString::number(step) + extension;

  QFile file(calcDir + QDir::separator() + fileName);
  if(file.open(IO_ReadOnly))
  {
    QTextStream stream(&file);
    //while(!stream.atEnd())
    //  result += stream.readLine();
    result = QStringList::split("\n", stream.read(), true);
  }
  return result;
}

std::vector<double> Calculation::loadFromPunch(const QString code, const unsigned int numValues, const unsigned int fieldSize, const unsigned int fieldsPerLine)
{
  std::vector<double> result;

  QFile file(calcDir + QDir::separator() + calcName + ".pun");
  if(!file.open(IO_ReadOnly))
    return result; // An empty result indicates failure

  qDebug("reading data from " + file.name());
  qDebug("looking for code |****" + code+ "|");

  QTextStream stream(&file);
  int codeLength = 4 + code.length();  
  while(stream.readLine().left(codeLength) != QString("****" + code) && !stream.atEnd())
    ;
  
  if(stream.atEnd())
  {
    result.clear();
    return result;
  }
  
  QString line;
  for(unsigned int i = 0; i < numValues/fieldsPerLine; i++)
  {
    line = stream.readLine(); // contains fieldsPerLine values
    for(unsigned int j = 0; j < 8; j++)
      result.push_back(line.mid(j*fieldSize,fieldSize).toDouble());

    if(stream.atEnd())
    {
      result.clear();
      return result;
    }
  }
  line = stream.readLine(); 
  for(unsigned int j = 0; j < numValues % fieldsPerLine; j++)
    result.push_back(line.mid(j*fieldSize,fieldSize).toDouble());

  assert(result.size() == numValues);
  return result;
}

---