BOSampleStepper.C 42.4 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 64 65 66 67
      cout << "<timing name=\""
           << setw(15) << (*i).first << "\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << 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 390 391 392 393 394 395
      if ( iter > 0 && ionic_stepper )
      {
        ionic_stepper->compute_v(energy,fion);
      }
      // at this point, positions r0, velocities v0 and forces fion are
      // consistent
      double ekin_ion = 0.0, temp_ion = 0.0;
      if ( ionic_stepper )
      {
        ekin_ion = ionic_stepper->ekin();
        temp_ion = ionic_stepper->temp();
      }
396

397 398
      // print positions, velocities and forces at time t0
      if ( onpe0 )
399
      {
400 401 402
        cout << "<atomset>" << endl;
        cout << atoms.cell();
        for ( int is = 0; is < atoms.atom_list.size(); is++ )
403
        {
404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420
          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
421
        }
422
        cout << "</atomset>" << endl;
Francois Gygi committed
423 424 425 426 427 428 429
        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
430
        cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
431
             << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
432
        cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
433
             << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
434
        cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
435
             << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
436
        cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
437
             << atoms.cell().volume() << " </unit_cell_volume>" << endl;
438 439 440 441 442 443

        // 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
444
        cout << setprecision(8);
445 446
        cout << "  <econst> " << enthalpy+ekin_ion+ekin_stepper
             << " </econst>\n";
447 448
        cout << "  <ekin_ion> " << ekin_ion << " </ekin_ion>\n";
        cout << "  <temp_ion> " << temp_ion << " </temp_ion>\n";
Francois Gygi committed
449
      }
450

451
      if ( atoms_move )
452
      {
453
        if ( s_.constraints.size() > 0 )
454
        {
455 456 457 458 459
          s_.constraints.compute_forces(ionic_stepper->r0(), fion);
          if ( onpe0 )
          {
            s_.constraints.list_constraints(cout);
          }
460
        }
461 462 463
        // move atoms to new position: r0 <- r0 + v0*dt + dt2/m * fion
        ionic_stepper->compute_r(energy,fion);
        ef_.atoms_moved();
464
      }
465

466
      if ( compute_stress )
467
      {
468 469
        compute_sigma();
        print_stress();
470

471 472
        if ( cell_moves )
        {
473
          cell_stepper->compute_new_cell(enthalpy,sigma,fion);
474 475 476

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

478 479 480
          ef_.cell_moved();
          ef_.atoms_moved(); // modifications of the cell also move ions
        }
481
      }
482
    } // if !gs_only
483

484
    // Recalculate ground state wavefunctions
485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506
#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
507
    // wavefunction extrapolation
Francois Gygi committed
508
    if ( atoms_move && extrapolate_wf )
509
    {
Francois Gygi committed
510
      for ( int ispin = 0; ispin < nspin; ispin++ )
511 512 513
      {
        for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
        {
514
          if ( ntc_extrapolation )
515
          {
516 517 518 519 520 521 522
            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 )
523
            {
524
              // copy c on cv
525
              for ( int i = 0; i < len; i++ )
526
              {
527 528 529 530 531
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
532
              }
533 534 535
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
536 537 538 539 540
            }
            else if ( iter == 1 )
            {
              s_.wfv->align(s_.wf);
              for ( int i = 0; i < len; i++ )
541
              {
542 543 544 545 546
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
                cmm[i] = xm;
547
              }
548 549 550 551 552 553 554 555 556 557 558 559
              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++ )
560
              {
561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581
                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];
582
              }
583 584 585 586 587 588 589 590 591 592 593 594 595 596 597
            }
            // 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++ )
              {
598 599 600 601 602
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
603
              }
604 605 606
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 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
            }
            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 )
            {
664
              // copy c to cv
665 666
              for ( int i = 0; i < len; i++ )
              {
667 668 669 670
                const double x = c[i];
                const double v = cv[i];
                c[i] = x + dt * v;
                cv[i] = x;
671
              }
672 673 674 675 676 677
              //tmap["lowdin"].start();
              //s_.wf.sd(ispin,ikp)->lowdin();
              //tmap["lowdin"].stop();
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
678 679 680
            }
            else
            {
681 682 683 684
              tmap["align"].start();
              s_.wfv->align(s_.wf);
              tmap["align"].stop();

685 686 687 688 689 690 691 692
              // 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;
              }
693 694 695
              //tmap["ortho_align"].start();
              //s_.wf.sd(ispin,ikp)->ortho_align(*s_.wfv->sd(ispin,ikp));
              //tmap["ortho_align"].stop();
696 697 698 699

              //tmap["riccati"].start();
              //s_.wf.sd(ispin,ikp)->riccati(*s_.wfv->sd(ispin,ikp));
              //tmap["riccati"].stop();
700 701 702 703

              tmap["lowdin"].start();
              s_.wf.sd(ispin,ikp)->lowdin();
              tmap["lowdin"].stop();
704
            }
705 706 707
          }
        }
      }
Francois Gygi committed
708
    } // atoms_move && extrapolate_wf
709

Francois Gygi committed
710
    // do nitscf self-consistent iterations, each with nite electronic steps
711
    if ( wf_stepper != 0 )
712
    {
713
      wf_stepper->preprocess();
Francois Gygi committed
714
      if ( anderson_charge_mixing )
715
        mixer.restart();
Francois Gygi committed
716

717
      double ehart, ehart_m;
718 719
      bool scf_converged = false;
      int itscf = 0;
720
      double etotal = 0.0, etotal_m = 0.0, etotal_mm = 0.0;
721

722
      while ( !scf_converged && itscf < nitscf_ )
723
      {
724
        if ( nite_ > 0 && onpe0 )
725
          cout << "  BOSampleStepper: start scf iteration" << endl;
726

727
        // compute new density in cd_.rhog
Francois Gygi committed
728
        tmap["charge"].start();
729 730 731 732
        if ( itscf==0 && initial_atomic_density )
          cd_.update_rhor();
        else
          cd_.update_density();
Francois Gygi committed
733
        tmap["charge"].stop();
734

Francois Gygi committed
735
        // charge mixing
736
        if ( nite_ > 0 )
737
        {
738
          if ( itscf == 0 )
Francois Gygi committed
739
          {
740
            // at first scf iteration, nothing to mix
741
            // memorize rhog in rhog_in for next iteration
742

743
            for ( int ispin = 0; ispin < nspin; ispin++ )
744 745
              for ( int i = 0; i < ng; i++ )
                rhog_in[i+ng*ispin] = cd_.rhog[ispin][i];
746 747 748
          }
          else
          {
749 750
            // itscf > 0
            // compute unscreened correction drhog
751
            for ( int ispin = 0; ispin < nspin; ispin++ )
752
            {
753
              for ( int i = 0; i < ng; i++ )
754
              {
755
                drhog[i+ng*ispin] = (cd_.rhog[ispin][i]-rhog_in[i+ng*ispin]);
756 757
              }
            }
758 759 760 761

            const double alpha = s_.ctrl.charge_mix_coeff;
            // Anderson acceleration
            if ( anderson_charge_mixing )
762
            {
763
              // row weighting of LS calculation
764
              for ( int ispin = 0; ispin < nspin; ispin++ )
765
              {
766 767
                for ( int i = 0; i < ng; i++ )
                  drhog[i+ng*ispin] /= wls[i];
768
              }
769

770 771 772 773
              const Context * const kpctxt = s_.wf.kpcontext();
              if ( kpctxt->mycol() == 0 )
              {
                // use AndersonMixer on first column only and bcast results
774 775 776 777 778
                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);
779 780 781
              }
              else
              {
782 783 784
                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);
785
              }
786

787
              for ( int ispin = 0; ispin < nspin; ispin++ )
788
              {
789 790 791 792 793
                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];
794
              }
795 796
            }
            else
797
            {
798
              for ( int ispin = 0; ispin < nspin; ispin++ )
799
              {
800 801
                for ( int i = 0; i < ng; i++ )
                  rhog_in[i+ng*ispin] += alpha * drhog[i+ng*ispin] * wkerker[i];
802
              }
803
            }
804

805
            for ( int ispin = 0; ispin < nspin; ispin++ )
806
            {
807 808
              for ( int i = 0; i < ng; i++ )
                cd_.rhog[ispin][i] = rhog_in[i+ng*ispin];
809
            }
810
            cd_.update_rhor();
811
          }
812
        } // if nite_ > 0
Francois Gygi committed
813

Francois Gygi committed
814
        tmap["update_vhxc"].start();
815
        ef_.update_vhxc(compute_stress);
Francois Gygi committed
816
        tmap["update_vhxc"].stop();
817

818
        // reset stepper only if multiple non-selfconsistent steps
819
        if ( nite_ > 0 ) wf_stepper->preprocess();
820

821
        // non-self-consistent loop
822 823
        // repeat until the change in eigenvalue_sum is smaller than a
        // fraction of the change in Hartree energy in the last scf iteration
824
        bool nonscf_converged = false;
825 826
        if ( itscf > 0 )
          ehart_m = ehart;
827
        ehart = ef_.ehart();
828 829 830
        double delta_ehart = 0.0;
        if ( itscf > 0 )
          delta_ehart = fabs(ehart - ehart_m);
831
        // if ( onpe0 && nite_ > 0 )
Francois Gygi committed
832
        //   cout << " delta_ehart = " << delta_ehart << endl;
833
        int ite = 0;
834
        double energy, etotal_int;
835 836 837 838 839 840

        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
841
        {
Francois Gygi committed
842
          tmap["energy"].start();
843
          energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
844
          tmap["energy"].stop();
845

Francois Gygi committed
846 847 848 849
          // 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
850 851
          // Note: since the hamiltonian is hermitian and dwf=H*wf
          // the dot product in the following line is real
852

853 854 855
          if ( ite > 0 )
            eigenvalue_sum_m = eigenvalue_sum;

856
          eigenvalue_sum = real(s_.wf.dot(dwf));
Francois Gygi committed
857
          if ( onpe0 )
858
            cout << "  <eigenvalue_sum>  "
Francois Gygi committed
859
                 << eigenvalue_sum << " </eigenvalue_sum>" << endl;
860

Francois Gygi committed
861
          tmap["wf_update"].start();
Francois Gygi committed
862
          wf_stepper->update(dwf);
Francois Gygi committed
863
          tmap["wf_update"].stop();
864

865
          // compare delta_eig_sum only after first iteration
866 867 868
          if ( ite > 0 )
          {
            double delta_eig_sum = fabs(eigenvalue_sum - eigenvalue_sum_m);
869 870
            nonscf_converged |= (delta_eig_sum < 0.01 * delta_ehart);
#ifdef DEBUG
871
            if ( onpe0 )
872 873 874 875 876 877
            {
              cout << " BOSampleStepper::step delta_eig_sum: "
                   << delta_eig_sum << endl;
              cout << " BOSampleStepper::step: delta_ehart: "
                   << delta_ehart << endl;
            }
878
#endif
879 880 881
          }
          ite++;
        }
882

Francois Gygi committed
883 884 885
        // subspace diagonalization
        if ( compute_eigvec || s_.ctrl.wf_diag == "EIGVAL" )
        {
Francois Gygi committed
886
          tmap["energy"].start();
887
          ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
888 889
          tmap["energy"].stop();
          tmap["wfdiag"].start();
Francois Gygi committed
890
          s_.wf.diag(dwf,compute_eigvec);
Francois Gygi committed
891
          tmap["wfdiag"].stop();
892
          if ( onpe0 )
893
          {
Francois Gygi committed
894
            cout << "<eigenset>" << endl;
895 896 897 898 899
            // print eigenvalues
            for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
            {
              for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
              {
900 901 902
                const int nst = wf.sd(ispin,ikp)->nst();
                const double eVolt = 2.0 * 13.6058;
                cout <<    "  <eigenvalues spin=\"" << ispin
903 904 905
                     << "\" kpoint=\""
                     << setprecision(8)
                     << wf.sd(ispin,ikp)->kpoint()
906 907 908
                     << "\" weight=\""
                     << setprecision(8)
                     << wf.weight(ikp)
909 910
                     << "\" n=\"" << nst << "\">" << endl;
                for ( int i = 0; i < nst; i++ )
911
                {
912 913 914
                  cout << setw(12) << setprecision(5)
                       << wf.sd(ispin,ikp)->eig(i)*eVolt;
                  if ( i%5 == 4 ) cout << endl;
915
                }
916 917
                if ( nst%5 != 0 ) cout << endl;
                cout << "  </eigenvalues>" << endl;
918
              }
919
            }
Francois Gygi committed
920
            cout << "</eigenset>" << endl;
921
          }
922
        }
923

924
        // update occupation numbers if fractionally occupied states
925
        // compute weighted sum of eigenvalues
926 927 928
        // default value if no fractional occupation
        double w_eigenvalue_sum = 2.0 * eigenvalue_sum;

929
        if ( fractional_occ )
Francois Gygi committed
930
        {
931
          wf.update_occ(s_.ctrl.fermi_temp);
932
#if 0
Francois Gygi committed
933
          if ( onpe0 )
Francois Gygi committed
934
          {
935 936 937
            cout << "  Wavefunction entropy: " << wf_entropy << endl;
            cout << "  Entropy contribution to free energy: "
                 << - wf_entropy * s_.ctrl.fermi_temp * boltz << endl;
Francois Gygi committed
938
          }
939
#endif
940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957
          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() +
958
                     ef_.esr() - ef_.eself() + ef_.dxc() + ef_.ets();
959 960 961 962 963 964 965 966 967 968
#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
969
        }
970 971 972 973 974 975 976 977 978 979
#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";
        }

980 981
        etotal_mm = etotal_m;
        etotal_m = etotal;
982
        etotal = etotal_int;
983

984
        if ( nite_ > 0 && onpe0 )
985
          cout << "  BOSampleStepper: end scf iteration" << endl;
986

987
        // delta_etotal = interval containing etotal, etotal_m and etotal_mm
988
        double delta_etotal = fabs(etotal - etotal_m);
989 990
        delta_etotal = max(delta_etotal,fabs(etotal - etotal_mm));
        delta_etotal = max(delta_etotal,fabs(etotal_m - etotal_mm));
991 992 993
        scf_converged |= (delta_etotal < s_.ctrl.scf_tol);
        itscf++;
      } // while scf
994

995 996
      if ( compute_mlwf || compute_mlwfc )
      {
Francois Gygi committed
997
        tmap["mlwf"].start();
998
        for ( int ispin = 0; ispin < nspin; ispin++ )
999
        {
1000
          mlwft[ispin]->update();
1001
          mlwft[ispin]->compute_transform();
1002
        }
1003 1004

        if ( compute_mlwf )
1005
        {
1006
          for ( int ispin = 0; ispin < nspin; ispin++ )
1007
          {
1008 1009
            SlaterDet& sd = *(wf.sd(ispin,0));
            mlwft[ispin]->apply_transform(sd);
1010 1011
          }
        }
1012 1013 1014

        if ( onpe0 )
        {
1015
          D3vector edipole_sum;
1016
          cout << "<mlwfs>" << endl;
1017
          for ( int ispin = 0; ispin < nspin; ispin++ )
1018
          {
1019
            SlaterDet& sd = *(wf.sd(ispin,0));
1020
            cout << " <mlwfset spin=\"" << ispin
1021
                 << "\" size=\"" << sd.nst() << "\">" << endl;
1022 1023 1024
            double total_spread[6];
            for ( int j = 0; j < 6; j++ )
               total_spread[j] = 0.0;
1025 1026
            for ( int i = 0; i < sd.nst(); i++ )
            {
1027
              D3vector ctr = mlwft[ispin]->center(i);
1028 1029 1030 1031 1032 1033 1034
              double spi[6];
              for (int j=0; j<3; j++)
              {
                spi[j] = mlwft[ispin]->spread2(i,j);
                total_spread[j] += spi[j];
              }

1035 1036 1037 1038 1039 1040
              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
1041 1042 1043 1044
                   << " \" spread=\" "
                   << setw(12) << spi[0]
                   << setw(12) << spi[1]
                   << setw(12) << spi[2] << " \"/>"
1045 1046 1047
                   << endl;
            }

1048 1049 1050 1051
            cout << " <total_spread> ";
            for ( int j = 0; j < 3; j++ )
              cout << setw(10) << total_spread[j];
            cout << " </total_spread>" << endl;
1052 1053
            D3vector edipole = mlwft[ispin]->dipole();
            cout << " <electronic_dipole spin=\"" << ispin << "\"> " << edipole
1054
                 << " </electronic_dipole>" << endl;
1055
            edipole_sum += edipole;
1056
            cout << " </mlwfset>" << endl;
1057
          }
1058 1059 1060 1061 1062 1063 1064
          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;
1065
          cout << "</mlwfs>" << endl;
1066
        }
Francois Gygi committed
1067
        tmap["mlwf"].stop();
1068 1069
      }

1070 1071 1072 1073
      // If GS calculation only, print energy and atomset at end of iterations
      if ( gs_only )
      {
        tmap["charge"].start();
1074
        cd_.update_density();
1075 1076
        tmap["charge"].stop();

Francois Gygi committed
1077
        tmap["update_vhxc"].start();
1078
        ef_.update_vhxc(compute_stress);
Francois Gygi committed
1079
        tmap["update_vhxc"].stop();
1080
        const bool compute_forces = true;
Francois Gygi committed
1081
        tmap["energy"].start();
1082
        ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
Francois Gygi committed
1083
        tmap["energy"].stop();
1084

Francois Gygi committed
1085
        if ( onpe0 )
1086
        {
Francois Gygi committed
1087
          cout << ef_;
1088 1089
          if ( ef_.el_enth() )
            cout << *ef_.el_enth();
1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112
          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
1113 1114 1115 1116 1117 1118 1119
          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
1120
          cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
1121
               << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
1122
          cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
1123
               << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
1124
          cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
1125
               << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
1126
          cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
1127
               << atoms.cell().volume() << " </unit_cell_volume>" << endl;
1128 1129 1130 1131 1132 1133 1134
          if ( compute_stress )
          {
            compute_sigma();
            print_stress();
          }
        }
      }
1135
      wf_stepper->postprocess();
1136
    }
1137 1138 1139 1140
    else
    {
      // wf_stepper == 0, wf_dyn == LOCKED
      // evaluate and print energy
Francois Gygi committed
1141
      tmap["charge"].start();
Francois Gygi committed
1142
      cd_.update_density();
Francois Gygi committed
1143
      tmap["charge"].stop();
Francois Gygi committed
1144
      tmap["update_vhxc"].start();
1145
      ef_.update_vhxc(compute_stress);
Francois Gygi committed
1146 1147
      tmap["update_vhxc"].stop();
      tmap["energy"].start();
Francois Gygi committed
1148
      ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
1149
      tmap["energy"].stop();
Francois Gygi committed
1150
      if ( onpe0 )
1151 1152
      {
        cout << ef_;
1153 1154 1155
        if ( ef_.el_enth() )
          cout << *ef_.el_enth();

1156 1157
      }
    }
1158

1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172
    if ( atoms_move )
      s_.constraints.update_constraints(dt);

    // execute commands in iter_cmd if defined
    if ( !s_.ctrl.iter_cmd.empty() )
    {
      if ( iter % s_.ctrl.iter_cmd_period == 0 )
      {
        // command must be terminated with \n
        istringstream cmdstream(s_.ctrl.iter_cmd + "\n");
        s_.ui->processCmds(cmdstream,"[iter_cmd]",true);
      }
    }

1173 1174 1175 1176 1177 1178
    // 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
1179
    if ( onpe0 )
1180
    {
1181 1182 1183 1184 1185
      cout << "  <timing name=\"iteration\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
      cout << "</iteration>" << endl;
1186 1187
    }
  } // for iter
1188

Francois Gygi committed
1189
  if ( atoms_move )
1190
  {
1191 1192
    // compute ionic forces at last position to update velocities
    // consistently with last position
Francois Gygi committed
1193
    tmap["charge"].start();
Francois Gygi committed
1194
    cd_.update_density();
Francois Gygi committed
1195
    tmap["charge"].stop();
1196

Francois Gygi committed
1197
    tmap["update_vhxc"].start();
1198
    ef_.update_vhxc(compute_stress);
Francois Gygi committed
1199
    tmap["update_vhxc"].stop();
Francois Gygi committed
1200
    const bool compute_forces = true;
Francois Gygi committed
1201
    tmap["energy"].start();
1202
    double energy =
1203
      ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
Francois Gygi committed
1204
    tmap["energy"].stop();
1205

Francois Gygi committed
1206 1207
    ionic_stepper->compute_v(energy,fion);
    // positions r0 and velocities v0 are consistent
1208
  }
1209

1210 1211 1212 1213 1214
  if ( atoms_move && extrapolate_wf )
  {
    // compute wavefunction velocity after last iteration
    // s_.wfv contains the previous wavefunction

Francois Gygi committed
1215
    tmap["align"].start();
1216
    s_.wfv->align(s_.wf);
Francois Gygi committed
1217
    tmap["align"].stop();
1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262

    for ( int ispin = 0; ispin < s_.wf.nspin(); ispin++ )
    {
      for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
      {
        double* c = (double*)