BOSampleStepper.C 42.7 KB
Newer Older
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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
15 16 17 18 19 20 21 22 23 24 25
// BOSampleStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "BOSampleStepper.h"
#include "EnergyFunctional.h"
#include "SlaterDet.h"
#include "Basis.h"
#include "WavefunctionStepper.h"
#include "SDWavefunctionStepper.h"
#include "PSDWavefunctionStepper.h"
26
#include "PSDAWavefunctionStepper.h"
27
#include "JDWavefunctionStepper.h"
28
#include "UserInterface.h"
29
#include "Preconditioner.h"
30
#include "SDIonicStepper.h"
Francois Gygi committed
31
#include "SDAIonicStepper.h"
Francois Gygi committed
32
#include "CGIonicStepper.h"
33
#include "MDIonicStepper.h"
Francois Gygi committed
34
#include "BMDIonicStepper.h"
35
#include "SDCellStepper.h"
36
#include "CGCellStepper.h"
37
#include "AndersonMixer.h"
38
#include "MLWFTransform.h"
39
#include "D3tensor.h"
40 41 42 43 44 45

#include <iostream>
#include <iomanip>
using namespace std;

////////////////////////////////////////////////////////////////////////////////
46 47
BOSampleStepper::BOSampleStepper(Sample& s, int nitscf, int nite) :
  SampleStepper(s), cd_(s.wf), ef_(s,cd_),
48 49
  dwf(s.wf), wfv(s.wfv), nitscf_(nitscf), nite_(nite),
  initial_atomic_density(false) {}
50

Francois Gygi committed
51 52 53 54 55 56 57 58 59 60 61 62
////////////////////////////////////////////////////////////////////////////////
BOSampleStepper::~BOSampleStepper()
{
  for ( TimerMap::iterator i = tmap.begin(); i != tmap.end(); i++ )
  {
    double time = (*i).second.real();
    double tmin = time;
    double tmax = time;
    s_.ctxt_.dmin(1,1,&tmin,1);
    s_.ctxt_.dmax(1,1,&tmax,1);
    if ( s_.ctxt_.myproc()==0 )
    {
63
      cout << "<timing name=\""
64 65 66
           << (*i).first << "\""
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
67
           << endl;
Francois Gygi committed
68 69 70 71
    }
  }
}

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
////////////////////////////////////////////////////////////////////////////////
void BOSampleStepper::initialize_density(void)
{
  // initialize cd_ with a sum of atomic densities

  double atom_radius[] =
  {
    0.00,  0.80,  0.59,  3.16,  2.12,  // null H  He  Li  Be
    1.64,  1.27,  1.06,  0.91,  0.79,  //  B   C   N   O   F
    0.72,  3.59,  2.74,  2.23,  2.10,  // Ne  Na  Mg  Al  Si
    1.85,  1.66,  1.49,  1.34,  4.59,  //  P   S  Cl  Ar   K
    3.67,  3.48,  3.33,  3.23,  3.14,  // Ca  Sc  Ti   V  Cr
    3.04,  2.95,  2.87,  2.82,  2.74,  // Mn  Fe  Co  Ni  Cu
    2.68,  2.57,  2.36,  2.15,  1.95,  // Zn  Ga  Ge  As  Se
    1.78,  1.66,  5.01,  4.14,  4.01,  // Br  Kr  Rb  Sr   Y
    3.89,  3.74,  3.59,  3.46,  3.36,  // Zr  Nb  Mo  Tc  Ru
    3.27,  3.19,  3.12,  3.04,  2.95,  // Rh  Pd  Ag  Cd  In
    2.74,  2.51,  2.32,  2.17,  2.04,  // Sn  Sb  Te   I  Xe
    5.63,  4.78,  0.00,  0.00,  4.67,  // Cs  Ba  La  Ce  Pr
    3.89,  3.87,  4.50,  4.37,  4.40,  // Nd  Pm  Sm  Eu  Gd
    4.25,  4.31,  0.00,  4.27,  4.20,  // Tb  Dy  Ho  Er  Tm
    4.20,  4.10,  3.93,  3.78,  3.65,  // Yb  Lu  Hf  Ta   W
    3.55,  3.50,  3.40,  3.34,  3.29,  // Re  Os  Ir  Pt  Au
    3.23,  2.95,  2.91,  2.70,  2.55,  // Hg  Tl  Pb  Bi  Po
    4.00,  2.27,  4.00,  4.00,  4.00,  // At  Rn  Fr  Ra  Ac
    4.00,  4.00,  4.00,  4.00,  4.00,  // Th  Pa   U  Np  Pu
    4.00,  4.00,  4.00,  4.00,  4.00,  // Am  Cm  Bk  Cf  Es
    4.00,  4.00,  4.00,  4.00          // Fm  Md  No  Lr
  };

  const AtomSet& atoms = s_.atoms;
  const Basis* const vbasis = cd_.vbasis();
  const int ngloc = vbasis->localsize();
  vector<vector<complex<double> > > rhops;
  const int nsp = atoms.nsp();
  rhops.resize(nsp);
  vector<complex<double> > rhopst(ngloc);
  const double * const g2 = vbasis->g2_ptr();

  for ( int is = 0; is < nsp; is++ )
  {
    rhops[is].resize(ngloc);
    Species *s = atoms.species_list[is];
    const int zval = s->zval();
    const int atomic_number = s->atomic_number();
    assert(atomic_number < sizeof(atom_radius)/sizeof(double));
    // scaling factor 2.0 in next line is empirically adjusted
    double rc = 2.0 * atom_radius[atomic_number];
    for ( int ig = 0; ig < ngloc; ig++ )
    {
      double arg = 0.25 * rc * rc * g2[ig];
      rhops[is][ig] = zval * exp( -arg );
    }
  }

  vector<vector<double> > tau0;
  tau0.resize(nsp);
  for ( int is = 0; is < nsp; is++ )
    tau0.resize(3*atoms.na(is));
  atoms.get_positions(tau0);
  StructureFactor sf;
  sf.init(tau0,*vbasis);
  sf.update(tau0,*vbasis);

  memset( (void*)&rhopst[0], 0, 2*ngloc*sizeof(double) );
  for ( int is = 0; is < nsp; is++ )
  {
    complex<double> *s = &sf.sfac[is][0];
    for ( int ig = 0; ig < ngloc; ig++ )
    {
      const complex<double> sg = s[ig];
      rhopst[ig] += sg * rhops[is][ig];
    }
  }

147
  // Initialize charge equally for both spins
148
  cd_.rhog[0] = rhopst;
149 150 151 152 153 154 155 156 157
  if ( cd_.rhog.size() == 2 )
  {
    assert(cd_.rhog[0].size()==cd_.rhog[1].size());
    for ( int i = 0; i < cd_.rhog[0].size(); i++ )
    {
      cd_.rhog[0][i] = 0.5 * rhopst[i];
      cd_.rhog[1][i] = 0.5 * rhopst[i];
    }
  }
158 159 160
  initial_atomic_density = true;
}

161
////////////////////////////////////////////////////////////////////////////////
162
void BOSampleStepper::step(int niter)
163
{
164 165
  const Context& ctxt = s_.ctxt_;
  const bool onpe0 = ctxt.onpe0();
Francois Gygi committed
166

167
  const bool anderson_charge_mixing = ( s_.ctrl.charge_mix_ndim > 0 );
Francois Gygi committed
168

169 170 171 172 173
  // determine whether eigenvectors must be computed
  // eigenvectors are computed if explicitly requested with wf_diag==T
  // or if the SlaterDet has fractionally occupied states
  const bool fractional_occ = (s_.wf.nel() != 2 * s_.wf.nst());
  const bool compute_eigvec = fractional_occ || s_.ctrl.wf_diag == "T";
174 175
  const bool compute_mlwf = s_.ctrl.wf_diag == "MLWF";
  const bool compute_mlwfc = s_.ctrl.wf_diag == "MLWFC";
176
  enum ortho_type { GRAM, LOWDIN, ORTHO_ALIGN, RICCATI };
177

178 179
  AtomSet& atoms = s_.atoms;
  Wavefunction& wf = s_.wf;
Francois Gygi committed
180
  const int nspin = wf.nspin();
181

182
  const double dt = s_.ctrl.dt;
183

184 185
  const string wf_dyn = s_.ctrl.wf_dyn;
  const string atoms_dyn = s_.ctrl.atoms_dyn;
186
  const string cell_dyn = s_.ctrl.cell_dyn;
187

188
  const bool extrapolate_wf = ( atoms_dyn == "MD" );
Francois Gygi committed
189

190
  const bool ntc_extrapolation =
191
    s_.ctrl.debug.find("NTC_EXTRAPOLATION") != string::npos;
192
  const bool asp_extrapolation =
193
    s_.ctrl.debug.find("ASP_EXTRAPOLATION") != string::npos;
194 195

  Wavefunction* wfmm;
Francois Gygi committed
196 197
  if ( extrapolate_wf && ( ntc_extrapolation || asp_extrapolation ) )
    wfmm = new Wavefunction(wf);
198

Francois Gygi committed
199 200
  // Next lines: special value of niter = 0: GS calculation only
  const bool atoms_move = ( niter > 0 && atoms_dyn != "LOCKED" );
201
  const bool compute_stress = ( s_.ctrl.stress == "ON" );
Francois Gygi committed
202 203
  const bool cell_moves = ( niter > 0 && compute_stress &&
                            cell_dyn != "LOCKED" );
204 205
  // GS-only calculation:
  const bool gs_only = !atoms_move && !cell_moves;
206

207
  Timer tm_iter;
208

209 210
  Preconditioner prec(wf,ef_,s_.ctrl.ecutprec);

211 212
  WavefunctionStepper* wf_stepper = 0;
  if ( wf_dyn == "SD" )
213 214 215
  {
    const double emass = s_.ctrl.emass;
    double dt2bye = (emass == 0.0) ? 0.5 / wf.ecut() : dt*dt/emass;
216

217 218 219 220 221 222 223 224
    // divide dt2bye by facs coefficient if stress == ON
    const double facs = 2.0;
    if ( s_.ctrl.stress == "ON" )
    {
      dt2bye /= facs;
    }
    wf_stepper = new SDWavefunctionStepper(wf,dt2bye,tmap);
  }
225
  else if ( wf_dyn == "PSD" )
226
    wf_stepper = new PSDWavefunctionStepper(wf,prec,tmap);
227
  else if ( wf_dyn == "PSDA" )
228
    wf_stepper = new PSDAWavefunctionStepper(wf,prec,tmap);
229
  else if ( wf_dyn == "JD" )
230
    wf_stepper = new JDWavefunctionStepper(wf,prec,ef_,tmap);
231

232
  // wf_stepper == 0 indicates that wf_dyn == LOCKED
233

234 235 236
  IonicStepper* ionic_stepper = 0;
  if ( atoms_dyn == "SD" )
    ionic_stepper = new SDIonicStepper(s_);
Francois Gygi committed
237 238
  else if ( atoms_dyn == "SDA" )
    ionic_stepper = new SDAIonicStepper(s_);
Francois Gygi committed
239 240
  else if ( atoms_dyn == "CG" )
    ionic_stepper = new CGIonicStepper(s_);
241 242
  else if ( atoms_dyn == "MD" )
    ionic_stepper = new MDIonicStepper(s_);
Francois Gygi committed
243 244
  else if ( atoms_dyn == "BMD" )
    ionic_stepper = new BMDIonicStepper(s_);
245

Francois Gygi committed
246 247
  if ( ionic_stepper )
    ionic_stepper->setup_constraints();
248

249 250 251
  CellStepper* cell_stepper = 0;
  if ( cell_dyn == "SD" )
    cell_stepper = new SDCellStepper(s_);
252 253
  else if ( cell_dyn == "CG" )
    cell_stepper = new CGCellStepper(s_);
254

255
  // Allocate wavefunction velocity if not available
Francois Gygi committed
256
  if ( atoms_move && extrapolate_wf )
257 258
  {
    if ( s_.wfv == 0 )
259
    {
260
      s_.wfv = new Wavefunction(wf);
261 262 263 264
      s_.wfv->clear();
    }
  }

265
  vector<MLWFTransform*> mlwft(nspin);
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280

  if ( compute_mlwf || compute_mlwfc )
  {
    // MLWF can be computed at the gamma point only
    // There must be a single k-point, and it must be gamma
    if ( wf.nkp() > 1 || ( wf.nkp()==1 && wf.kpoint(0) != D3vector(0,0,0) ) )
    {
      if ( onpe0 )
      {
        cout << " BOSampleStepper::step: MLWF can be computed at k=0 only"
             << endl;
        cout << " BOSampleStepper::step: cannot run" << endl;
      }
      return;
    }
281 282 283

    for ( int ispin = 0; ispin < nspin; ispin++ )
      mlwft[ispin] = new MLWFTransform(*wf.sd(ispin,0));
284
  }
285

286 287 288 289 290
  // Charge mixing variables: include both spins in the same vector
  if ( nspin > 1 ) assert(cd_.rhog[0].size()==cd_.rhog[1].size());
  const int ng = cd_.rhog[0].size();
  vector<complex<double> > rhog_in(nspin*ng), drhog(nspin*ng),
    rhobar(nspin*ng), drhobar(nspin*ng);
291
  const int anderson_ndim = s_.ctrl.charge_mix_ndim;
292
  vector<double> wkerker(ng), wls(ng);
293

294 295 296 297
  // Anderson charge mixer: include both spins in the same vector
  // Factor of 2: complex coeffs stored as double
  MPI_Comm vcomm = cd_.vcomm();
  AndersonMixer mixer(2*nspin*ng,anderson_ndim,&vcomm);
298 299 300 301 302

  // compute Kerker preconditioning
  // real space Kerker cutoff in a.u.
  const double rc_Kerker = s_.ctrl.charge_mix_rcut;
  const double *const g2 = cd_.vbasis()->g2_ptr();
303 304 305 306 307 308 309 310 311 312 313

  // define q1 cutoff for row weighting of LS charge mixing
  // Use rc1 = 3 a.u. default cutoff
  double rc1 = 3.0;
  // check if override from the debug variable
  // use: set debug RC1 <value>
  if ( s_.ctrl.debug.find("RC1") != string::npos )
  {
    istringstream is(s_.ctrl.debug);
    string s;
    is >> s >> rc1;
314 315
    if ( onpe0 )
      cout << " override rc1 value: rc1 = " << rc1 << endl;
316 317 318 319 320 321
    assert(rc1 >= 0.0);
  }

  if ( rc1 != 0.0 )
  {
    const double q1 = 2.0 * M_PI / rc1;
322
    for ( int i = 0; i < wls.size(); i++ )
323 324 325 326 327 328 329 330 331
    {
      if ( g2[i] != 0.0 )
        wls[i] = sqrt(g2[i] / ( g2[i] + q1*q1 ));
      else
        wls[i] = 1.0;
    }
  }
  else
  {
332
    for ( int i = 0; i < wls.size(); i++ )
333 334 335
      wls[i] = 1.0;
  }

336 337 338 339
  if ( rc_Kerker > 0.0 )
  {
    const double q0_kerker = 2 * M_PI / rc_Kerker;
    const double q0_kerker2 = q0_kerker * q0_kerker;
340
    for ( int i = 0; i < wkerker.size(); i++ )
341
      wkerker[i] = g2[i] / ( g2[i] + q0_kerker2 );
342 343 344
  }
  else
  {
345
    for ( int i = 0; i < wkerker.size(); i++ )
346 347 348
      wkerker[i] = 1.0;
  }

Francois Gygi committed
349 350
  // Next line: special case of niter=0: compute GS only
  for ( int iter = 0; iter < max(niter,1); iter++ )
351 352
  {
    // ionic iteration
353

354
    tm_iter.start();
355

Francois Gygi committed
356
    if ( onpe0 )
357
      cout << "<iteration count=\"" << iter+1 << "\">\n";
358

Francois Gygi committed
359
    // compute energy and ionic forces using existing wavefunction
Francois Gygi committed
360

361 362 363
    if ( !gs_only )
    {
      tmap["charge"].start();
364
      cd_.update_density();
365
      tmap["charge"].stop();
366

Francois Gygi committed
367
      tmap["update_vhxc"].start();
368
      ef_.update_vhxc(compute_stress);
Francois Gygi committed
369
      tmap["update_vhxc"].stop();
370
      const bool compute_forces = true;
Francois Gygi committed
371
      tmap["energy"].start();
372 373
      double energy =
        ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
Francois Gygi committed
374
      tmap["energy"].stop();
375
      double enthalpy = ef_.enthalpy();
376

377
      if ( onpe0 )
378
      {
379
        cout << ef_;
380 381
        if ( ef_.el_enth() )
          cout << *ef_.el_enth();
382
      }
383

384 385 386 387 388 389
      if ( iter > 0 && ionic_stepper )
      {
        ionic_stepper->compute_v(energy,fion);
      }
      // at this point, positions r0, velocities v0 and forces fion are
      // consistent
390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413

      // execute commands in iter_cmd if defined
      if ( !iter_cmd_.empty() )
      {
        if ( iter % iter_cmd_period_ == 0 )
        {
          // copy positions and velocities from IonicStepper to AtomSet
          if ( ionic_stepper )
          {
            ionic_stepper->set_positions();
            ionic_stepper->set_velocities();
          }
          // command must be terminated with \n
          istringstream cmdstream(iter_cmd_ + "\n");
          s_.ui->processCmds(cmdstream,"[iter_cmd]",true);
          // copy positions and velocities back from AtomSet
          if ( ionic_stepper )
          {
            ionic_stepper->get_positions();
            ionic_stepper->get_velocities();
          }
        }
      }

414 415 416 417 418 419
      double ekin_ion = 0.0, temp_ion = 0.0;
      if ( ionic_stepper )
      {
        ekin_ion = ionic_stepper->ekin();
        temp_ion = ionic_stepper->temp();
      }
420

421 422
      // print positions, velocities and forces at time t0
      if ( onpe0 )
423
      {
424 425 426
        cout << "<atomset>" << endl;
        cout << atoms.cell();
        for ( int is = 0; is < atoms.atom_list.size(); is++ )
427
        {
428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444
          int i = 0;
          for ( int ia = 0; ia < atoms.atom_list[is].size(); ia++ )
          {
            Atom* pa = atoms.atom_list[is][ia];
            cout << "  <atom name=\"" << pa->name() << "\""
                 << " species=\"" << pa->species()
                 << "\">\n"
                 << "    <position> " << pa->position() << " </position>\n"
                 << "    <velocity> " << pa->velocity() << " </velocity>\n"
                 << "    <force> "
                 << fion[is][i] << " "
                 << fion[is][i+1] << " "
                 << fion[is][i+2]
                 << " </force>\n";
            cout << "  </atom>" << endl;
            i += 3;
          }
Francois Gygi committed
445
        }
446
        cout << "</atomset>" << endl;
Francois Gygi committed
447 448 449 450 451 452 453
        cout << setprecision(6);
        cout << "<unit_cell_a_norm> " << atoms.cell().a_norm(0)
             << " </unit_cell_a_norm>" << endl;
        cout << "<unit_cell_b_norm> " << atoms.cell().a_norm(1)
             << " </unit_cell_b_norm>" << endl;
        cout << "<unit_cell_c_norm> " << atoms.cell().a_norm(2)
             << " </unit_cell_c_norm>" << endl;
Francois Gygi committed
454
        cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
455
             << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
456
        cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
457
             << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
458
        cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
459
             << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
460
        cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
461
             << atoms.cell().volume() << " </unit_cell_volume>" << endl;
462 463 464 465 466 467

        // include the kinetic energy of the stepper
        // e.g. to include thermostat contributions
        double ekin_stepper;
        if ( ionic_stepper != 0 )
          ekin_stepper = ionic_stepper->ekin_stepper();
Francois Gygi committed
468
        cout << setprecision(8);
469 470
        cout << "  <econst> " << enthalpy+ekin_ion+ekin_stepper
             << " </econst>\n";
471 472
        cout << "  <ekin_ion> " << ekin_ion << " </ekin_ion>\n";
        cout << "  <temp_ion> " << temp_ion << " </temp_ion>\n";
Francois Gygi committed
473
      }
474

475
      if ( atoms_move )
476
      {
477
        if ( s_.constraints.size() > 0 )
478
        {
479 480 481 482 483
          s_.constraints.compute_forces(ionic_stepper->r0(), fion);
          if ( onpe0 )
          {
            s_.constraints.list_constraints(cout);
          }
484
        }
485 486 487
        // move atoms to new position: r0 <- r0 + v0*dt + dt2/m * fion
        ionic_stepper->compute_r(energy,fion);
        ef_.atoms_moved();
488
      }
489

490
      if ( compute_stress )
491
      {
492 493
        compute_sigma();
        print_stress();
494

495 496
        if ( cell_moves )
        {
497
          cell_stepper->compute_new_cell(enthalpy,sigma,fion);
498 499 500

          // Update cell
          cell_stepper->update_cell();
501

502 503 504
          ef_.cell_moved();
          ef_.atoms_moved(); // modifications of the cell also move ions
        }
505
      }
506
    } // if !gs_only
507

508
    // Recalculate ground state wavefunctions
509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530
#ifdef DEBUG
    for ( int ispin = 0; ispin < nspin; ispin++ )
    {
      for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
      {
        double sum = s_.wf.sd(ispin,ikp)->empty_row_error() ;
        if ( onpe0 )
        {
          cout.setf(ios::scientific,ios::floatfield);
          cout << " sd empty row error: ispin="
               << ispin << " ikp=" << ikp << "  " << sum << endl;
        }
        sum = s_.wf.sd(ispin,ikp)->g0_imag_error() ;
        if ( onpe0 )
        {
          cout.setf(ios::scientific,ios::floatfield);
          cout << " sd g0 imag error: ispin="
               << ispin << " ikp=" << ikp << "  " << sum << endl;
        }
      }
    }
#endif
531
    // wavefunction extrapolation
Francois Gygi committed
532
    if ( atoms_move && extrapolate_wf )
533
    {
Francois Gygi committed
534
      for ( int ispin = 0; ispin < nspin; ispin++ )
535 536 537
      {
        for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
        {
538
          if ( ntc_extrapolation )
539
          {
540 541 542 543 544 545 546
            double* c = (double*) s_.wf.sd(ispin,ikp)->c().cvalptr();
            double* cv = (double*) s_.wfv->sd(ispin,ikp)->c().cvalptr();
            double* cmm = (double*) wfmm->sd(ispin,ikp)->c().cvalptr();
            const int mloc = s_.wf.sd(ispin,ikp)->c().mloc();
            const int nloc = s_.wf.sd(ispin,ikp)->c().nloc();
            const int len = 2*mloc*nloc;
            if ( iter == 0 )
547
            {
548
              // copy c on cv
549
              for ( int i = 0; i < len; i++ )
550
              {
551 552 553 554 555
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
556
              }
557 558 559
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
560 561 562 563 564
            }
            else if ( iter == 1 )
            {
              s_.wfv->align(s_.wf);
              for ( int i = 0; i < len; i++ )
565
              {
566 567 568 569 570
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
                cmm[i] = xm;
571
              }
572 573 574 575 576 577 578 579 580 581 582 583
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
            }
            else
            {
              // align wf with wfmm before extrapolation
              // s_.wf.align(*wfmm);
              wfmm->align(s_.wf);

              // extrapolate
              for ( int i = 0; i < len; i++ )
584
              {
585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605
                const double x = c[i];   // current wf (scf converged) at t
                const double xm = cv[i]; // extrapolated wf at t
                const double xmm = cmm[i]; // extrapolated wf at t-dt
                c[i] = 2.0 * x - xmm;
                // save extrapolated value at t in cmm
                cmm[i] = xm;
              }

              // orthogonalize the extrapolated value
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
              //tmap["lowdin"].start();
              //s_.wf.sd(ispin,ikp)->lowdin();
              //tmap["lowdin"].stop();

              // c[i] now contains the extrapolated value
              // save a copy in cv[i]
              for ( int i = 0; i < len; i++ )
              {
                cv[i] = c[i];
606
              }
607 608 609 610 611 612 613 614 615 616 617 618 619 620 621
            }
            // c[i] is now ready for electronic iterations
          }
          else if ( asp_extrapolation )
          {
            double* c = (double*) s_.wf.sd(ispin,ikp)->c().cvalptr();
            double* cv = (double*) s_.wfv->sd(ispin,ikp)->c().cvalptr();
            double* cmm = (double*) wfmm->sd(ispin,ikp)->c().cvalptr();
            const int mloc = s_.wf.sd(ispin,ikp)->c().mloc();
            const int nloc = s_.wf.sd(ispin,ikp)->c().nloc();
            const int len = 2*mloc*nloc;
            if ( iter == 0 )
            {
              for ( int i = 0; i < len; i++ )
              {
622 623 624 625 626
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
627
              }
628 629 630
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687
            }
            else if ( iter == 1 )
            {
              //s_.wfv->align(s_.wf);
              for ( int i = 0; i < len; i++ )
              {
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
                cmm[i] = xm;
              }
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
            }
            else
            {
              // align wf with wfmm before extrapolation
              // s_.wf.align(*wfmm);
              // wfmm->align(s_.wf);

              // extrapolate
              for ( int i = 0; i < len; i++ )
              {
                const double x = c[i];   // current wf (scf converged) at t
                const double xm = cv[i]; // extrapolated wf at t
                const double xmm = cmm[i]; // extrapolated wf at t-dt
                const double asp_a1 = 0.5;
                c[i] = 2.0 * x - xm +
                       asp_a1 * ( x - 2.0 * xm + xmm );
                //c[i] = 2.5 * x - 2.0 * xm + 0.5 * xmm;
                cmm[i] = xm;
                cv[i] = x;
              }

              // orthogonalize the extrapolated value
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
              //tmap["lowdin"].start();
              //s_.wf.sd(ispin,ikp)->lowdin();
              //tmap["lowdin"].stop();

              // c[i] now contains the extrapolated value
            }
            // c[i] is now ready for electronic iterations
          }
          else // normal extrapolation
          {
            double* c = (double*) s_.wf.sd(ispin,ikp)->c().cvalptr();
            double* cv = (double*) s_.wfv->sd(ispin,ikp)->c().cvalptr();
            const int mloc = s_.wf.sd(ispin,ikp)->c().mloc();
            const int nloc = s_.wf.sd(ispin,ikp)->c().nloc();
            const int len = 2*mloc*nloc;
            if ( iter == 0 )
            {
688
              // copy c to cv
689 690
              for ( int i = 0; i < len; i++ )
              {
691 692 693 694
                const double x = c[i];
                const double v = cv[i];
                c[i] = x + dt * v;
                cv[i] = x;
695
              }
696 697 698 699 700 701
              //tmap["lowdin"].start();
              //s_.wf.sd(ispin,ikp)->lowdin();
              //tmap["lowdin"].stop();
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
702 703 704
            }
            else
            {
705 706 707 708
              tmap["align"].start();
              s_.wfv->align(s_.wf);
              tmap["align"].stop();

709 710 711 712 713 714 715 716
              // linear extrapolation
              for ( int i = 0; i < len; i++ )
              {
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
              }
717 718 719
              //tmap["ortho_align"].start();
              //s_.wf.sd(ispin,ikp)->ortho_align(*s_.wfv->sd(ispin,ikp));
              //tmap["ortho_align"].stop();
720 721 722 723

              //tmap["riccati"].start();
              //s_.wf.sd(ispin,ikp)->riccati(*s_.wfv->sd(ispin,ikp));
              //tmap["riccati"].stop();
724 725 726 727

              tmap["lowdin"].start();
              s_.wf.sd(ispin,ikp)->lowdin();
              tmap["lowdin"].stop();
728
            }
729 730 731
          }
        }
      }
Francois Gygi committed
732
    } // atoms_move && extrapolate_wf
733

Francois Gygi committed
734
    // do nitscf self-consistent iterations, each with nite electronic steps
735
    if ( wf_stepper != 0 )
736
    {
737
      wf_stepper->preprocess();
Francois Gygi committed
738
      if ( anderson_charge_mixing )
739
        mixer.restart();
Francois Gygi committed
740

741
      double ehart, ehart_m;
742 743
      bool scf_converged = false;
      int itscf = 0;
744
      double etotal = 0.0, etotal_m = 0.0, etotal_mm = 0.0;
745

746
      while ( !scf_converged && itscf < nitscf_ )
747
      {
748
        if ( nite_ > 0 && onpe0 )
749
          cout << "  BOSampleStepper: start scf iteration" << endl;
750

751
        // compute new density in cd_.rhog
Francois Gygi committed
752
        tmap["charge"].start();
753 754 755 756
        if ( itscf==0 && initial_atomic_density )
          cd_.update_rhor();
        else
          cd_.update_density();
Francois Gygi committed
757
        tmap["charge"].stop();
758

Francois Gygi committed
759
        // charge mixing
760
        if ( nite_ > 0 )
761
        {
762
          if ( itscf == 0 )
Francois Gygi committed
763
          {
764
            // at first scf iteration, nothing to mix
765
            // memorize rhog in rhog_in for next iteration
766

767
            for ( int ispin = 0; ispin < nspin; ispin++ )
768 769
              for ( int i = 0; i < ng; i++ )
                rhog_in[i+ng*ispin] = cd_.rhog[ispin][i];
770 771 772
          }
          else
          {
773 774
            // itscf > 0
            // compute unscreened correction drhog
775
            for ( int ispin = 0; ispin < nspin; ispin++ )
776
            {
777
              for ( int i = 0; i < ng; i++ )
778
              {
779
                drhog[i+ng*ispin] = (cd_.rhog[ispin][i]-rhog_in[i+ng*ispin]);
780 781
              }
            }
782 783 784 785

            const double alpha = s_.ctrl.charge_mix_coeff;
            // Anderson acceleration
            if ( anderson_charge_mixing )
786
            {
787
              // row weighting of LS calculation
788
              for ( int ispin = 0; ispin < nspin; ispin++ )
789
              {
790 791
                for ( int i = 0; i < ng; i++ )
                  drhog[i+ng*ispin] /= wls[i];
792
              }
793

794 795 796 797
              const Context * const kpctxt = s_.wf.kpcontext();
              if ( kpctxt->mycol() == 0 )
              {
                // use AndersonMixer on first column only and bcast results
798 799 800 801 802
                mixer.update((double*)&rhog_in[0], (double*)&drhog[0],
                             (double*)&rhobar[0], (double*)&drhobar[0]);
                const int n = 2*nspin*ng;
                kpctxt->dbcast_send('r',n,1,(double*)&rhobar[0],n);
                kpctxt->dbcast_send('r',n,1,(double*)&drhobar[0],n);
803 804 805
              }
              else
              {
806 807 808
                const int n = 2*nspin*ng;
                kpctxt->dbcast_recv('r',n,1,(double*)&rhobar[0],n,-1,0);
                kpctxt->dbcast_recv('r',n,1,(double*)&drhobar[0],n,-1,0);
809
              }
810

811
              for ( int ispin = 0; ispin < nspin; ispin++ )
812
              {
813 814 815 816 817
                for ( int i = 0; i < ng; i++ )
                  drhobar[i+ng*ispin] *= wls[i];
                for ( int i = 0; i < ng; i++ )
                  rhog_in[i+ng*ispin] = rhobar[i+ng*ispin] +
                    alpha * drhobar[i+ng*ispin] * wkerker[i];
818
              }
819 820
            }
            else
821
            {
822
              for ( int ispin = 0; ispin < nspin; ispin++ )
823
              {
824 825
                for ( int i = 0; i < ng; i++ )
                  rhog_in[i+ng*ispin] += alpha * drhog[i+ng*ispin] * wkerker[i];
826
              }
827
            }
828

829
            for ( int ispin = 0; ispin < nspin; ispin++ )
830
            {
831 832
              for ( int i = 0; i < ng; i++ )
                cd_.rhog[ispin][i] = rhog_in[i+ng*ispin];
833
            }
834
            cd_.update_rhor();
835
          }
836
        } // if nite_ > 0
Francois Gygi committed
837

Francois Gygi committed
838
        tmap["update_vhxc"].start();
839
        ef_.update_vhxc(compute_stress);
Francois Gygi committed
840
        tmap["update_vhxc"].stop();
841

842
        // reset stepper only if multiple non-selfconsistent steps
843
        if ( nite_ > 0 ) wf_stepper->preprocess();
844

845
        // non-self-consistent loop
846 847
        // repeat until the change in eigenvalue_sum is smaller than a
        // fraction of the change in Hartree energy in the last scf iteration
848
        bool nonscf_converged = false;
849 850
        if ( itscf > 0 )
          ehart_m = ehart;
851
        ehart = ef_.ehart();
852 853 854
        double delta_ehart = 0.0;
        if ( itscf > 0 )
          delta_ehart = fabs(ehart - ehart_m);
855
        // if ( onpe0 && nite_ > 0 )
Francois Gygi committed
856
        //   cout << " delta_ehart = " << delta_ehart << endl;
857
        int ite = 0;
858
        double energy, etotal_int;
859 860 861 862 863 864

        double eigenvalue_sum, eigenvalue_sum_m = 0.0;
        // if nite == 0: do 1 iteration, no screening in charge mixing
        // if nite > 0: do nite iterations, use screening in charge mixing
        //
        while ( !nonscf_converged && ite < max(nite_,1) )
Francois Gygi committed
865
        {
Francois Gygi committed
866
          tmap["energy"].start();
867
          energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
868
          tmap["energy"].stop();
869

Francois Gygi committed
870 871 872 873
          // compute the sum of eigenvalues (with fixed weight)
          // to measure convergence of the subspace update
          // compute trace of the Hamiltonian matrix Y^T H Y
          // scalar product of Y and (HY): tr Y^T (HY) = sum_ij Y_ij (HY)_ij
874 875
          // Note: since the hamiltonian is hermitian and dwf=H*wf
          // the dot product in the following line is real
876

877 878 879
          if ( ite > 0 )
            eigenvalue_sum_m = eigenvalue_sum;

880
          eigenvalue_sum = real(s_.wf.dot(dwf));
Francois Gygi committed
881
          if ( onpe0 )
882
            cout << "  <eigenvalue_sum>  "
Francois Gygi committed
883
                 << eigenvalue_sum << " </eigenvalue_sum>" << endl;
884

Francois Gygi committed
885
          tmap["wf_update"].start();
Francois Gygi committed
886
          wf_stepper->update(dwf);
Francois Gygi committed
887
          tmap["wf_update"].stop();
888

889
          // compare delta_eig_sum only after first iteration
890 891 892
          if ( ite > 0 )
          {
            double delta_eig_sum = fabs(eigenvalue_sum - eigenvalue_sum_m);
893 894
            nonscf_converged |= (delta_eig_sum < 0.01 * delta_ehart);
#ifdef DEBUG
895
            if ( onpe0 )
896 897 898 899 900 901
            {
              cout << " BOSampleStepper::step delta_eig_sum: "
                   << delta_eig_sum << endl;
              cout << " BOSampleStepper::step: delta_ehart: "
                   << delta_ehart << endl;
            }
902
#endif
903 904 905
          }
          ite++;
        }
906

Francois Gygi committed
907 908 909
        // subspace diagonalization
        if ( compute_eigvec || s_.ctrl.wf_diag == "EIGVAL" )
        {
Francois Gygi committed
910
          tmap["energy"].start();
911
          ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
912 913
          tmap["energy"].stop();
          tmap["wfdiag"].start();
Francois Gygi committed
914
          s_.wf.diag(dwf,compute_eigvec);
Francois Gygi committed
915
          tmap["wfdiag"].stop();
916
          if ( onpe0 )
917
          {
Francois Gygi committed
918
            cout << "<eigenset>" << endl;
919 920 921 922 923
            // print eigenvalues
            for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
            {
              for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
              {
924 925 926
                const int nst = wf.sd(ispin,ikp)->nst();
                const double eVolt = 2.0 * 13.6058;
                cout <<    "  <eigenvalues spin=\"" << ispin
927 928 929
                     << "\" kpoint=\""
                     << setprecision(8)
                     << wf.sd(ispin,ikp)->kpoint()
930 931 932
                     << "\" weight=\""
                     << setprecision(8)
                     << wf.weight(ikp)
933 934
                     << "\" n=\"" << nst << "\">" << endl;
                for ( int i = 0; i < nst; i++ )
935
                {
936 937 938
                  cout << setw(12) << setprecision(5)
                       << wf.sd(ispin,ikp)->eig(i)*eVolt;
                  if ( i%5 == 4 ) cout << endl;
939
                }
940 941
                if ( nst%5 != 0 ) cout << endl;
                cout << "  </eigenvalues>" << endl;
942
              }
943
            }
Francois Gygi committed
944
            cout << "</eigenset>" << endl;
945
          }
946
        }
947

948
        // update occupation numbers if fractionally occupied states
949
        // compute weighted sum of eigenvalues
950 951 952
        // default value if no fractional occupation
        double w_eigenvalue_sum = 2.0 * eigenvalue_sum;

953
        if ( fractional_occ )
Francois Gygi committed
954
        {
955
          wf.update_occ(s_.ctrl.fermi_temp);
956
#if 0
Francois Gygi committed
957
          if ( onpe0 )
Francois Gygi committed
958
          {
959 960 961
            cout << "  Wavefunction entropy: " << wf_entropy << endl;
            cout << "  Entropy contribution to free energy: "
                 << - wf_entropy * s_.ctrl.fermi_temp * boltz << endl;
Francois Gygi committed
962
          }
963
#endif
964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981
          w_eigenvalue_sum = 0.0;
          for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
          {
            for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
            {
              const int nst = wf.sd(ispin,ikp)->nst();
              const double wkp = wf.weight(ikp);
              for ( int n = 0; n < nst; n++ )
              {
                const double occ = wf.sd(ispin,ikp)->occ(n);
                w_eigenvalue_sum += wkp * occ * wf.sd(ispin,ikp)->eig(n);
              }
            }
          }
        }

        // Harris-Foulkes estimate of the total energy
        etotal_int = w_eigenvalue_sum - ef_.ehart_e() + ef_.ehart_p() +
982
                     ef_.esr() - ef_.eself() + ef_.dxc() + ef_.ets();
983 984 985 986 987 988 989 990 991 992
#ifdef DEBUG
        if ( onpe0 )
        {
          cout << setprecision(8);
          cout << "w_eigenvalue_sum = " << setw(15) << w_eigenvalue_sum << endl;
          cout << "ef.dxc()         = " << setw(15) << ef_.dxc() << endl;
          cout << "ef.ehart()       = " << setw(15) << ef_.ehart() << endl;
          cout << "ef.ehart_e()     = " << setw(15) << ef_.ehart_e() << endl;
          cout << "ef.ehart_ep()    = " << setw(15) << ef_.ehart_ep() << endl;
          cout << "ef.esr()         = " << setw(15) << ef_.esr() << endl;
Francois Gygi committed
993
        }
994 995 996 997 998 999 1000 1001 1002 1003
#endif

        if ( onpe0 )
        {
          cout.setf(ios::fixed,ios::floatfield);
          cout.setf(ios::right,ios::adjustfield);
          cout << "  <etotal_int>  " << setprecision(8) << setw(15)
               << etotal_int << " </etotal_int>\n";
        }

1004 1005
        etotal_mm = etotal_m;
        etotal_m = etotal;
1006
        etotal = etotal_int;
1007

1008
        if ( nite_ > 0 && onpe0 )
1009
          cout << "  BOSampleStepper: end scf iteration" << endl;
1010

1011
        // delta_etotal = interval containing etotal, etotal_m and etotal_mm
1012
        double delta_etotal = fabs(etotal - etotal_m);
1013 1014
        delta_etotal = max(delta_etotal,fabs(etotal - etotal_mm));
        delta_etotal = max(delta_etotal,fabs(etotal_m - etotal_mm));
1015 1016 1017
        scf_converged |= (delta_etotal < s_.ctrl.scf_tol);
        itscf++;
      } // while scf
1018

1019 1020
      if ( compute_mlwf || compute_mlwfc )
      {
Francois Gygi committed
1021
        tmap["mlwf"].start();
1022
        for ( int ispin = 0; ispin < nspin; ispin++ )
1023
        {
1024
          mlwft[ispin]->update();
1025
          mlwft[ispin]->compute_transform();
1026
        }
1027 1028

        if ( compute_mlwf )
1029
        {
1030
          for ( int ispin = 0; ispin < nspin; ispin++ )
1031
          {
1032 1033
            SlaterDet& sd = *(wf.sd(ispin,0));
            mlwft[ispin]->apply_transform(sd);
1034 1035
          }
        }
1036 1037 1038

        if ( onpe0 )
        {
1039
          D3vector edipole_sum;
1040
          cout << "<mlwfs>" << endl;
1041
          for ( int ispin = 0; ispin < nspin; ispin++ )
1042
          {
1043
            SlaterDet& sd = *(wf.sd(ispin,0));
1044
            cout << " <mlwfset spin=\"" << ispin
1045
                 << "\" size=\"" << sd.nst() << "\">" << endl;
1046 1047 1048
            double total_spread[6];
            for ( int j = 0; j < 6; j++ )
               total_spread[j] = 0.0;
1049 1050
            for ( int i = 0; i < sd.nst(); i++ )
            {
1051
              D3vector ctr = mlwft[ispin]->center(i);
1052 1053 1054 1055 1056 1057 1058
              double spi[6];
              for (int j=0; j<3; j++)
              {
                spi[j] = mlwft[ispin]->spread2(i,j);
                total_spread[j] += spi[j];
              }

1059 1060 1061 1062 1063 1064
              cout.setf(ios::fixed, ios::floatfield);
              cout.setf(ios::right, ios::adjustfield);
              cout << "   <mlwf center=\"" << setprecision(6)
                   << setw(12) << ctr.x
                   << setw(12) << ctr.y
                   << setw(12) << ctr.z
1065 1066 1067 1068
                   << " \" spread=\" "
                   << setw(12) << spi[0]
                   << setw(12) << spi[1]
                   << setw(12) << spi[2] << " \"/>"
1069 1070 1071
                   << endl;
            }

1072 1073 1074 1075
            cout << " <total_spread> ";
            for ( int j = 0; j < 3; j++ )
              cout << setw(10) << total_spread[j];
            cout << " </total_spread>" << endl;
1076 1077
            D3vector edipole = mlwft[ispin]->dipole();
            cout << " <electronic_dipole spin=\"" << ispin << "\"> " << edipole
1078
                 << " </electronic_dipole>" << endl;
1079
            edipole_sum += edipole;
1080
            cout << " </mlwfset>" << endl;
1081
          }
1082 1083 1084 1085 1086 1087 1088
          D3vector idipole = atoms.dipole();
          cout << " <ionic_dipole> " << idipole
               << " </ionic_dipole>" << endl;
          cout << " <total_dipole> " << idipole + edipole_sum
               << " </total_dipole>" << endl;
          cout << " <total_dipole_length> " << length(idipole + edipole_sum)
               << " </total_dipole_length>" << endl;
1089
          cout << "</mlwfs>" << endl;
1090
        }
Francois Gygi committed
1091
        tmap["mlwf"].stop();
1092 1093
      }

1094 1095 1096 1097
      // If GS calculation only, print energy and atomset at end of iterations
      if ( gs_only )
      {
        tmap["charge"].start();
1098
        cd_.update_density();
1099 1100
        tmap["charge"].stop();

Francois Gygi committed
1101
        tmap["update_vhxc"].start();
1102
        ef_.update_vhxc(compute_stress);
Francois Gygi committed
1103
        tmap["update_vhxc"].stop();
1104
        const bool compute_forces = true;
Francois Gygi committed
1105
        tmap["energy"].start();
1106
        ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
Francois Gygi committed
1107
        tmap["energy"].stop();
1108

Francois Gygi committed
1109
        if ( onpe0 )
1110
        {
Francois Gygi committed
1111
          cout << ef_;
1112 1113
          if ( ef_.el_enth() )
            cout << *ef_.el_enth();
1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136
          cout << "<atomset>" << endl;
          cout << atoms.cell();
          for ( int is = 0; is < atoms.atom_list.size(); is++ )
          {
            int i = 0;
            for ( int ia = 0; ia < atoms.atom_list[is].size(); ia++ )
            {
              Atom* pa = atoms.atom_list[is][ia];
              cout << "  <atom name=\"" << pa->name() << "\""
                   << " species=\"" << pa->species()
                   << "\">\n"
                   << "    <position> " << pa->position() << " </position>\n"
                   << "    <velocity> " << pa->velocity() << " </velocity>\n"
                   << "    <force> "
                   << fion[is][i] << " "
                   << fion[is][i+1] << " "
                   << fion[is][i+2]
                   << " </force>\n";
              cout << "  </atom>" << endl;
              i += 3;
            }
          }
          cout << "</atomset>" << endl;
Francois Gygi committed
1137 1138 1139 1140 1141 1142 1143
          cout << setprecision(6);
          cout << "<unit_cell_a_norm> " << atoms.cell().a_norm(0)
               << " </unit_cell_a_norm>" << endl;
          cout << "<unit_cell_b_norm> " << atoms.cell().a_norm(1)
               << " </unit_cell_b_norm>" << endl;
          cout << "<unit_cell_c_norm> " << atoms.cell().a_norm(2)
               << " </unit_cell_c_norm>" << endl;
Francois Gygi committed
1144
          cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
1145
               << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
1146
          cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
1147
               << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
1148
          cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
1149
               << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
1150
          cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
1151
               << atoms.cell().volume() << " </unit_cell_volume>" << endl;
1152 1153 1154 1155 1156 1157 1158
          if ( compute_stress )
          {
            compute_sigma();
            print_stress();
          }
        }
      }
1159
      wf_stepper->postprocess();
1160
    }
1161 1162 1163 1164
    else
    {
      // wf_stepper == 0, wf_dyn == LOCKED
      // evaluate and print energy
Francois Gygi committed
1165
      tmap["charge"].start();
Francois Gygi committed
1166
      cd_.update_density();
Francois Gygi committed
1167
      tmap["charge"].stop();
Francois Gygi committed
1168
      tmap["update_vhxc"].start();
1169
      ef_.update_vhxc(compute_stress);
Francois Gygi committed
1170 1171
      tmap["update_vhxc"].stop();
      tmap["energy"].start();
Francois Gygi committed
1172
      ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
1173
      tmap["energy"].stop();
Francois Gygi committed
1174
      if ( onpe0 )
1175 1176
      {
        cout << ef_;
1177 1178 1179
        if ( ef_.el_enth() )
          cout << *ef_.el_enth();

1180 1181
      }
    }
1182

1183 1184 1185
    if ( atoms_move )
      s_.constraints.update_constraints(dt);

1186 1187 1188 1189 1190 1191
    // print iteration time
    double time = tm_iter.real();
    double tmin = time;
    double tmax = time;
    s_.ctxt_.dmin(1,1,&tmin,1);
    s_.ctxt_.dmax(1,1,&tmax,1);
Francois Gygi committed
1192
    if ( onpe0 )
1193
    {
1194
      cout << "  <timing name=\"iteration\""
1195 1196
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
1197 1198
           << endl;
      cout << "</iteration>" << endl;
1199 1200
    }
  } // for iter
1201

Francois Gygi committed
1202
  if ( atoms_move )
1203
  {
1204 1205
    // compute ionic forces at last position to update velocities
    // consistently with last position
Francois Gygi committed
1206
    tmap["charge"].start();
Francois Gygi committed
1207
    cd_.update_density();
Francois Gygi committed
1208
    tmap["charge"].stop();
1209

Francois Gygi committed
1210
    tmap["update_vhxc"].start();
1211
    ef_.update_vhxc(compute_stress);
Francois Gygi committed
1212
    tmap["update_vhxc"].stop();
Francois Gygi committed
1213
    const bool compute_forces = true;
Francois Gygi committed
1214
    tmap["energy"].start();
1215
    double energy =
1216
      ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
Francois Gygi committed
1217
    tmap["energy"].stop();
1218

Francois Gygi committed
1219 1220
    ionic_stepper->compute_v(energy,fion);
    // positions r0 and velocities v0 are consistent
1221
  }
1222

1223 1224 1225 1226 1227
  if ( atoms_move && extrapolate_wf )
  {
    // compute wavefunction velocity after last iteration
    // s_.wfv contains the previous wavefunction

Francois Gygi committed
1228
    tmap["align"].start();
1229
    s_.wfv->align(s_.wf);
Francois Gygi committed
1230
    tmap