RunCmd.C 2.56 KB
Newer Older
Francois Gygi committed
1 2
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
3 4 5 6
// Copyright (c) 2008 The Regents of the University of California
//
// This file is part of Qbox
//
Francois Gygi committed
7 8
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
Francois Gygi committed
9 10 11 12 13 14
// the License, or (at your option) any later version.
// See the file COPYING in the root directory of this distribution
// or <http://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////////////
//
Francois Gygi committed
15 16 17 18 19 20 21
// RunCmd.C:
//
////////////////////////////////////////////////////////////////////////////////

#include "RunCmd.h"
#include<iostream>
using namespace std;
22 23
#include "BOSampleStepper.h"
#include "CPSampleStepper.h"
Francois Gygi committed
24 25 26 27 28 29

#include<ctime>
#include<cassert>

int RunCmd::action(int argc, char **argv)
{
30
  if ( argc < 2 || argc > 5)
Francois Gygi committed
31 32
  {
    if ( ui->onpe0() )
33
    {
34 35 36
      cout << " use: run [-atomic_density] niter" << endl;
      cout << "      run [-atomic_density] niter nitscf" << endl;
      cout << "      run [-atomic_density] niter nitscf nite" << endl;
37
    }
Francois Gygi committed
38 39
    return 1;
  }
40

Francois Gygi committed
41 42 43
  if ( s->wf.nst() == 0 )
  {
    if ( ui->onpe0() )
44
      cout << " RunCmd: no states, cannot run" << endl;
Francois Gygi committed
45 46
    return 1;
  }
47 48 49
  if ( s->wf.ecut() == 0.0 )
  {
    if ( ui->onpe0() )
50
      cout << " RunCmd: ecut = 0.0, cannot run" << endl;
51 52
    return 1;
  }
53 54 55 56 57 58
  if ( s->wf.cell().volume() == 0.0 )
  {
    if ( ui->onpe0() )
      cout << " RunCmd: volume = 0.0, cannot run" << endl;
    return 1;
  }
59

60
  SampleStepper* stepper;
61

62 63 64 65 66 67 68 69 70 71
  int iarg = 1;
  bool atomic_density = false;
  if ( !strcmp(argv[iarg],"-atomic_density") )
  {
    atomic_density = true;
    iarg++;
    argc--;
  }

  int niter = atoi(argv[iarg]);
72
  int nite = 0;
Francois Gygi committed
73
  int nitscf = 1;
74 75
  if ( argc == 3 )
  {
Francois Gygi committed
76
    // run niter nitscf
77
    nitscf = atoi(argv[iarg+1]);
Francois Gygi committed
78 79 80 81
  }
  else if ( argc == 4 )
  {
    // run niter nitscf nite
82 83
    nitscf = atoi(argv[iarg+1]);
    nite = atoi(argv[iarg+2]);
84
  }
Francois Gygi committed
85 86 87

  s->extforces.setup(s->atoms);

88
  if ( s->ctrl.wf_dyn == "MD" )
Francois Gygi committed
89
    stepper = new CPSampleStepper(*s);
90
  else
Francois Gygi committed
91
    stepper = new BOSampleStepper(*s,nitscf,nite);
92

93
  assert(stepper!=0);
94 95
  stepper->set_iter_cmd(s->ctrl.iter_cmd);
  stepper->set_iter_cmd_period(s->ctrl.iter_cmd_period);
96

97 98 99
  if ( atomic_density )
    stepper->initialize_density();

100
  s->wf.info(cout,"wavefunction");
101
  stepper->step(niter);
102

103 104 105 106 107 108 109 110
  // Delete wave function velocity if not using atoms_dyn = MD
  if ( s->ctrl.atoms_dyn != "MD" )
  {
    if ( s->wfv != 0 )
      delete s->wfv;
    s->wfv = 0;
  }

111
  delete stepper;
112

Francois Gygi committed
113 114
  return 0;
}