BOSampleStepper.C 44.9 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
  dwf(s.wf), wfv(s.wfv), nitscf_(nitscf), nite_(nite),
49
  update_density_first_(true), update_vh_(true), update_vxc_(true) {}
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
      string s = "name=\"" + (*i).first + "\"";
      cout << "<timing " << left << setw(22) << s
65 66
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
67
           << endl;
Francois Gygi committed
68 69 70 71
    }
  }
}

72 73 74 75 76 77 78
////////////////////////////////////////////////////////////////////////////////
void BOSampleStepper::initialize_density(void)
{
  // initialize cd_ with a sum of atomic densities

  double atom_radius[] =
  {
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
    0.000, 1.542, 0.931, 1.727, 1.636, // null H  He  Li  Be
    1.845, 1.538, 1.311, 1.141, 1.014, //  B   C   N   O   F
    0.983, 1.238, 1.347, 1.376, 2.151, // Ne  Na  Mg  Al  Si
    1.927, 1.733, 1.563, 1.429, 1.546, //  P   S  Cl  Ar   K
    1.663, 1.604, 1.516, 1.434, 1.377, // Ca  Sc  Ti   V  Cr
    1.310, 1.245, 1.201, 1.162, 1.098, // Mn  Fe  Co  Ni  Cu
    1.077, 1.331, 1.415, 2.015, 1.880, // Zn  Ga  Ge  As  Se
    1.749, 1.630, 1.705, 1.819, 1.794, // Br  Kr  Rb  Sr   Y
    1.728, 1.664, 1.589, 1.523, 1.461, // Zr  Nb  Mo  Tc  Ru
    1.410, 1.348, 1.306, 1.303, 1.554, // Rh  Pd  Ag  Cd  In
    1.609, 1.611, 1.530, 1.514, 1.464, // Sn  Sb  Te   I  Xe
    1.946, 1.967, 1.943, 1.930, 1.920, // Cs  Ba  La  Ce  Pr
    1.910, 1.900, 1.890, 1.880, 1.870, // Nd  Pm  Sm  Eu  Gd
    1.860, 1.850, 1.840, 1.830, 1.820, // Tb  Dy  Ho  Er  Tm
    1.810, 1.800, 1.701, 1.658, 1.606, // Yb  Lu  Hf  Ta   W
    1.550, 1.500, 1.446, 1.398, 1.355, // Re  Os  Ir  Pt  Au
    1.314, 1.624, 1.659, 1.634, 1.620, // Hg  Tl  Pb  Bi  Po
    1.600, 1.500, 1.600, 1.600, 1.600, // At  Rn  Fr  Ra  Ac
    1.600, 1.600, 1.600, 1.600, 1.600, // Th  Pa   U  Np  Pu
    1.600, 1.600, 1.600, 1.600, 1.600, // Am  Cm  Bk  Cf  Es
    1.600, 1.600, 1.600, 1.600         // Fm  Md  No  Lr
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
  };

  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));
118
    double rc = atom_radius[atomic_number];
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
    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];
    }
  }

146 147 148
  // Adjust G=0 component of the charge if net_charge is non-zero
  rhopst[0] += s_.wf.nel() - atoms.nel();

149
  // Initialize charge equally for both spins
150
  cd_.rhog[0] = rhopst;
151 152 153 154 155 156 157 158 159
  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];
    }
  }
160 161
  cd_.update_rhor();
  update_density_first_ = false;
He Ma committed
162 163
}

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

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

172 173 174 175 176
  // 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";
177 178
  const bool compute_mlwf = s_.ctrl.wf_diag == "MLWF";
  const bool compute_mlwfc = s_.ctrl.wf_diag == "MLWFC";
179
  enum ortho_type { GRAM, LOWDIN, ORTHO_ALIGN, RICCATI };
180

181 182
  const double gpa = 29421.5;

183 184
  AtomSet& atoms = s_.atoms;
  Wavefunction& wf = s_.wf;
Francois Gygi committed
185
  const int nspin = wf.nspin();
186

187
  const double dt = s_.ctrl.dt;
188

189 190
  const string wf_dyn = s_.ctrl.wf_dyn;
  const string atoms_dyn = s_.ctrl.atoms_dyn;
191
  const string cell_dyn = s_.ctrl.cell_dyn;
192

193
  const bool extrapolate_wf = ( atoms_dyn == "MD" );
Francois Gygi committed
194

195
  const bool ntc_extrapolation =
196
    s_.ctrl.debug.find("NTC_EXTRAPOLATION") != string::npos;
197
  const bool asp_extrapolation =
198
    s_.ctrl.debug.find("ASP_EXTRAPOLATION") != string::npos;
199 200

  Wavefunction* wfmm;
Francois Gygi committed
201 202
  if ( extrapolate_wf && ( ntc_extrapolation || asp_extrapolation ) )
    wfmm = new Wavefunction(wf);
203

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

212 213 214
  const double force_tol = s_.ctrl.force_tol;
  const double stress_tol = s_.ctrl.stress_tol;

215
  Timer tm_iter;
216

217 218
  Preconditioner prec(wf,ef_,s_.ctrl.ecutprec);

219 220
  WavefunctionStepper* wf_stepper = 0;
  if ( wf_dyn == "SD" )
221 222 223
  {
    const double emass = s_.ctrl.emass;
    double dt2bye = (emass == 0.0) ? 0.5 / wf.ecut() : dt*dt/emass;
224

225 226 227 228 229 230 231 232
    // 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);
  }
233
  else if ( wf_dyn == "PSD" )
234
    wf_stepper = new PSDWavefunctionStepper(wf,prec,tmap);
235
  else if ( wf_dyn == "PSDA" )
236
    wf_stepper = new PSDAWavefunctionStepper(wf,prec,tmap);
237
  else if ( wf_dyn == "JD" )
238
    wf_stepper = new JDWavefunctionStepper(wf,prec,ef_,tmap);
239

240
  // wf_stepper == 0 indicates that wf_dyn == LOCKED
241

242 243 244
  IonicStepper* ionic_stepper = 0;
  if ( atoms_dyn == "SD" )
    ionic_stepper = new SDIonicStepper(s_);
Francois Gygi committed
245 246
  else if ( atoms_dyn == "SDA" )
    ionic_stepper = new SDAIonicStepper(s_);
Francois Gygi committed
247 248
  else if ( atoms_dyn == "CG" )
    ionic_stepper = new CGIonicStepper(s_);
249 250
  else if ( atoms_dyn == "MD" )
    ionic_stepper = new MDIonicStepper(s_);
Francois Gygi committed
251 252
  else if ( atoms_dyn == "BMD" )
    ionic_stepper = new BMDIonicStepper(s_);
253

Francois Gygi committed
254 255
  if ( ionic_stepper )
    ionic_stepper->setup_constraints();
256

257 258 259
  CellStepper* cell_stepper = 0;
  if ( cell_dyn == "SD" )
    cell_stepper = new SDCellStepper(s_);
260 261
  else if ( cell_dyn == "CG" )
    cell_stepper = new CGCellStepper(s_);
262

263
  // Allocate wavefunction velocity if not available
Francois Gygi committed
264
  if ( atoms_move && extrapolate_wf )
265 266
  {
    if ( s_.wfv == 0 )
267
    {
268
      s_.wfv = new Wavefunction(wf);
269 270 271 272
      s_.wfv->clear();
    }
  }

273
  vector<MLWFTransform*> mlwft(nspin);
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288

  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;
    }
289 290 291

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

294 295 296 297 298
  // 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);
299
  const int anderson_ndim = s_.ctrl.charge_mix_ndim;
300
  vector<double> wkerker(ng), wls(ng);
301

302 303 304 305
  // 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);
306 307 308 309 310

  // 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();
311 312 313 314 315 316 317 318 319 320 321

  // 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;
322 323
    if ( onpe0 )
      cout << " override rc1 value: rc1 = " << rc1 << endl;
324 325 326 327 328 329
    assert(rc1 >= 0.0);
  }

  if ( rc1 != 0.0 )
  {
    const double q1 = 2.0 * M_PI / rc1;
330
    for ( int i = 0; i < wls.size(); i++ )
331 332 333 334 335 336 337 338 339
    {
      if ( g2[i] != 0.0 )
        wls[i] = sqrt(g2[i] / ( g2[i] + q1*q1 ));
      else
        wls[i] = 1.0;
    }
  }
  else
  {
340
    for ( int i = 0; i < wls.size(); i++ )
341 342 343
      wls[i] = 1.0;
  }

344 345 346 347
  if ( rc_Kerker > 0.0 )
  {
    const double q0_kerker = 2 * M_PI / rc_Kerker;
    const double q0_kerker2 = q0_kerker * q0_kerker;
348
    for ( int i = 0; i < wkerker.size(); i++ )
349
      wkerker[i] = g2[i] / ( g2[i] + q0_kerker2 );
350 351 352
  }
  else
  {
353
    for ( int i = 0; i < wkerker.size(); i++ )
354 355 356
      wkerker[i] = 1.0;
  }

357 358 359
  if ( onpe0 )
    cout << "<net_charge> " << atoms.nel()-wf.nel() << " </net_charge>\n";

Francois Gygi committed
360
  // Next line: special case of niter=0: compute GS only
361 362
  bool iter_done = false;
  for ( int iter = 0; iter < max(niter,1) && !iter_done; iter++ )
363 364
  {
    // ionic iteration
365

366
    tm_iter.start();
367

Francois Gygi committed
368
    if ( onpe0 )
369
      cout << "<iteration count=\"" << iter+1 << "\">\n";
370

Francois Gygi committed
371
    // compute energy and ionic forces using existing wavefunction
372 373
    double maxforce = 0.0;
    double maxstress = 0.0;
Francois Gygi committed
374

375 376 377
    if ( !gs_only )
    {
      tmap["charge"].start();
378
      cd_.update_density();
379
      tmap["charge"].stop();
380

Francois Gygi committed
381
      tmap["update_vhxc"].start();
382
      ef_.update_vhxc(compute_stress);
Francois Gygi committed
383
      tmap["update_vhxc"].stop();
384
      const bool compute_forces = true;
Francois Gygi committed
385
      tmap["energy"].start();
386 387
      double energy =
        ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
Francois Gygi committed
388
      tmap["energy"].stop();
389
      double enthalpy = ef_.enthalpy();
390

391 392 393 394 395 396 397 398 399 400
      if ( force_tol > 0.0 )
      {
        maxforce = 0.0;
        for ( int is = 0; is < fion.size(); is++ )
          for ( int i = 0; i < fion[is].size(); i++ )
            maxforce = max(maxforce, fabs(fion[is][i]));
      }

      if ( stress_tol > 0.0 )
      {
401
        compute_sigma();
402
        for ( int i = 0; i < sigma.size(); i++ )
403
          maxstress = max(maxstress, gpa*fabs(sigma[i]));
404 405
      }

406
      if ( onpe0 )
407
      {
408
        cout << cd_;
409
        cout << ef_;
410 411
        if ( ef_.el_enth() )
          cout << *ef_.el_enth();
412
      }
413

414 415 416 417 418 419
      if ( iter > 0 && ionic_stepper )
      {
        ionic_stepper->compute_v(energy,fion);
      }
      // at this point, positions r0, velocities v0 and forces fion are
      // consistent
420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443

      // 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();
          }
        }
      }

444 445 446 447 448 449
      double ekin_ion = 0.0, temp_ion = 0.0;
      if ( ionic_stepper )
      {
        ekin_ion = ionic_stepper->ekin();
        temp_ion = ionic_stepper->temp();
      }
450

451 452
      // print positions, velocities and forces at time t0
      if ( onpe0 )
453
      {
454 455 456
        cout << "<atomset>" << endl;
        cout << atoms.cell();
        for ( int is = 0; is < atoms.atom_list.size(); is++ )
457
        {
458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474
          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
475
        }
476
        cout << "</atomset>" << endl;
Francois Gygi committed
477 478 479 480 481 482 483
        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
484
        cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
485
             << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
486
        cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
487
             << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
488
        cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
489
             << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
490
        cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
491
             << atoms.cell().volume() << " </unit_cell_volume>" << endl;
492 493 494 495 496 497

        // 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
498
        cout << setprecision(8);
499 500
        cout << "  <econst> " << enthalpy+ekin_ion+ekin_stepper
             << " </econst>\n";
501 502
        cout << "  <ekin_ion> " << ekin_ion << " </ekin_ion>\n";
        cout << "  <temp_ion> " << temp_ion << " </temp_ion>\n";
Francois Gygi committed
503
      }
504

505
      if ( atoms_move )
506
      {
507
        if ( s_.constraints.size() > 0 )
508
        {
509 510 511 512 513
          s_.constraints.compute_forces(ionic_stepper->r0(), fion);
          if ( onpe0 )
          {
            s_.constraints.list_constraints(cout);
          }
514
        }
515 516 517
        // move atoms to new position: r0 <- r0 + v0*dt + dt2/m * fion
        ionic_stepper->compute_r(energy,fion);
        ef_.atoms_moved();
518
      }
519

520
      if ( compute_stress )
521
      {
522 523
        compute_sigma();
        print_stress();
524

525 526
        if ( cell_moves )
        {
527
          cell_stepper->compute_new_cell(enthalpy,sigma,fion);
528 529 530

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

532 533 534
          ef_.cell_moved();
          ef_.atoms_moved(); // modifications of the cell also move ions
        }
535
      }
536
    } // if !gs_only
537

538
    // Recalculate ground state wavefunctions
539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560
#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
561
    // wavefunction extrapolation
Francois Gygi committed
562
    if ( atoms_move && extrapolate_wf )
563
    {
Francois Gygi committed
564
      for ( int ispin = 0; ispin < nspin; ispin++ )
565 566 567
      {
        for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
        {
568
          if ( ntc_extrapolation )
569
          {
570 571 572 573 574 575 576
            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 )
577
            {
578
              // copy c on cv
579
              for ( int i = 0; i < len; i++ )
580
              {
581 582 583 584 585
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
586
              }
587 588 589
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
590 591 592
            }
            else if ( iter == 1 )
            {
593
              s_.wfv->sd(ispin,ikp)->align(*s_.wf.sd(ispin,ikp));
594
              for ( int i = 0; i < len; i++ )
595
              {
596 597 598 599 600
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
                cmm[i] = xm;
601
              }
602 603 604 605 606 607 608
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
            }
            else
            {
              // align wf with wfmm before extrapolation
609
              s_.wf.sd(ispin,ikp)->align(*wfmm->sd(ispin,ikp));
610 611 612

              // extrapolate
              for ( int i = 0; i < len; i++ )
613
              {
614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634
                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];
635
              }
636 637 638 639 640 641 642 643 644 645 646 647 648 649 650
            }
            // 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++ )
              {
651 652 653 654 655
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
656
              }
657 658 659
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
660 661 662
            }
            else if ( iter == 1 )
            {
663
              s_.wfv->sd(ispin,ikp)->align(*s_.wf.sd(ispin,ikp));
664 665 666 667 668 669 670 671 672 673 674 675 676 677 678
              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
679
              s_.wf.sd(ispin,ikp)->align(*wfmm->sd(ispin,ikp));
680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715

              // 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 )
            {
716
              // copy c to cv
717 718
              for ( int i = 0; i < len; i++ )
              {
719 720 721 722
                const double x = c[i];
                const double v = cv[i];
                c[i] = x + dt * v;
                cv[i] = x;
723
              }
724 725 726 727 728 729
              //tmap["lowdin"].start();
              //s_.wf.sd(ispin,ikp)->lowdin();
              //tmap["lowdin"].stop();
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
730 731 732
            }
            else
            {
733
              tmap["align"].start();
734
              s_.wfv->sd(ispin,ikp)->align(*s_.wf.sd(ispin,ikp));
735 736
              tmap["align"].stop();

737 738 739 740 741 742 743 744
              // 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;
              }
745 746 747
              //tmap["ortho_align"].start();
              //s_.wf.sd(ispin,ikp)->ortho_align(*s_.wfv->sd(ispin,ikp));
              //tmap["ortho_align"].stop();
748 749 750 751

              //tmap["riccati"].start();
              //s_.wf.sd(ispin,ikp)->riccati(*s_.wfv->sd(ispin,ikp));
              //tmap["riccati"].stop();
752 753 754 755

              tmap["lowdin"].start();
              s_.wf.sd(ispin,ikp)->lowdin();
              tmap["lowdin"].stop();
756
            }
757 758 759
          }
        }
      }
Francois Gygi committed
760
    } // atoms_move && extrapolate_wf
761

Francois Gygi committed
762
    // do nitscf self-consistent iterations, each with nite electronic steps
763
    if ( wf_stepper != 0 )
764
    {
765
      wf_stepper->preprocess();
Francois Gygi committed
766
      if ( anderson_charge_mixing )
767
        mixer.restart();
Francois Gygi committed
768

769
      double ehart, ehart_m;
770
      double delta_ehart;
771 772
      bool scf_converged = false;
      int itscf = 0;
773
      double etotal = 0.0, etotal_m = 0.0, etotal_mm = 0.0;
774

775
      while ( !scf_converged && itscf < nitscf_ )
776
      {
777
        if ( nite_ > 0 && onpe0 )
778
          cout << "  BOSampleStepper: start scf iteration" << endl;
779

780
        // update charge density
Francois Gygi committed
781
        tmap["charge"].start();
782 783 784 785
        // The density is updated at the first step if update_density_first_
        // is true.
        // It is always updated after the first step
        if ( ( update_density_first_ || itscf>0 ) )
786
          cd_.update_density();
Francois Gygi committed
787
        tmap["charge"].stop();
788

789 790 791
        if ( onpe0 )
          cout << cd_;

Francois Gygi committed
792
        // charge mixing
793
        if ( nite_ > 0 )
794
        {
795
          if ( itscf == 0 )
Francois Gygi committed
796
          {
797
            // at first scf iteration, nothing to mix
798
            // memorize rhog in rhog_in for next iteration
799

800
            for ( int ispin = 0; ispin < nspin; ispin++ )
801 802
              for ( int i = 0; i < ng; i++ )
                rhog_in[i+ng*ispin] = cd_.rhog[ispin][i];
803 804 805
          }
          else
          {
806 807
            // itscf > 0
            // compute unscreened correction drhog
808
            for ( int ispin = 0; ispin < nspin; ispin++ )
809
            {
810
              for ( int i = 0; i < ng; i++ )
811
              {
812
                drhog[i+ng*ispin] = (cd_.rhog[ispin][i]-rhog_in[i+ng*ispin]);
813 814
              }
            }
815 816 817 818

            const double alpha = s_.ctrl.charge_mix_coeff;
            // Anderson acceleration
            if ( anderson_charge_mixing )
819
            {
820
              // row weighting of LS calculation
821
              for ( int ispin = 0; ispin < nspin; ispin++ )
822
              {
823 824
                for ( int i = 0; i < ng; i++ )
                  drhog[i+ng*ispin] /= wls[i];
825
              }
826

827 828 829 830
              const Context * const kpctxt = s_.wf.kpcontext();
              if ( kpctxt->mycol() == 0 )
              {
                // use AndersonMixer on first column only and bcast results
831 832 833 834 835
                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);
836 837 838
              }
              else
              {
839 840 841
                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);
842
              }
843

844
              for ( int ispin = 0; ispin < nspin; ispin++ )
845
              {
846 847 848 849 850
                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];
851
              }
852 853
            }
            else
854
            {
855
              for ( int ispin = 0; ispin < nspin; ispin++ )
856
              {
857 858
                for ( int i = 0; i < ng; i++ )
                  rhog_in[i+ng*ispin] += alpha * drhog[i+ng*ispin] * wkerker[i];
859
              }
860
            }
861

862
            for ( int ispin = 0; ispin < nspin; ispin++ )
863
            {
864 865
              for ( int i = 0; i < ng; i++ )
                cd_.rhog[ispin][i] = rhog_in[i+ng*ispin];
866
            }
867
            cd_.update_rhor();
868
          }
869
        } // if nite_ > 0
Francois Gygi committed
870

871 872 873 874
        // update vhxc:
        // at first scf step:
        // - update both vh and vxc
        // at later steps:
875
        // - update depending of values of update_vh_ and update_vxc_
Francois Gygi committed
876
        tmap["update_vhxc"].start();
877 878
        if ( itscf == 0 )
          ef_.update_vhxc(compute_stress);
He Ma committed
879
        else
880
          ef_.update_vhxc(compute_stress, update_vh_, update_vxc_);
Francois Gygi committed
881
        tmap["update_vhxc"].stop();
882

883
        // reset stepper only if multiple non-selfconsistent steps
884
        if ( nite_ > 0 ) wf_stepper->preprocess();
885

886
        // non-self-consistent loop
887
        // repeat until the change in eigenvalue_sum is smaller than a
888 889
        // fraction delta_ratio of the change in Hartree energy delta_ehart
        // in the last scf iteration
890
        bool nonscf_converged = false;
891 892 893 894 895 896 897
        if ( itscf == 0 )
        {
          ehart = ef_.ehart();
          delta_ehart = 0.0;
        }
        else if ( itscf == 1 )
        {
898
          ehart_m = ehart;
899
          ehart = ef_.ehart();
900
          delta_ehart = fabs(ehart - ehart_m);
901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917
        }
        else
        {
          // itscf > 1
          // only allow decrease in delta_ehart
          ehart_m = ehart;
          ehart = ef_.ehart();
          delta_ehart = min(delta_ehart,fabs(ehart - ehart_m));
        }
#if DEBUG
        if ( onpe0 )
        {
          cout << " BOSampleStepper::step: delta_ehart: "
               << delta_ehart << endl;
        }
#endif

918
        int ite = 0;
Francois Gygi committed
919
        double etotal_int;
920 921 922 923 924 925

        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
926
        {
Francois Gygi committed
927
          tmap["energy"].start();
Francois Gygi committed
928
          ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
929
          tmap["energy"].stop();
930

Francois Gygi committed
931 932 933 934
          // 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
935 936
          // Note: since the hamiltonian is hermitian and dwf=H*wf
          // the dot product in the following line is real
937

938 939 940
          if ( ite > 0 )
            eigenvalue_sum_m = eigenvalue_sum;

941
          eigenvalue_sum = real(s_.wf.dot(dwf));
Francois Gygi committed
942
          if ( onpe0 )
943
            cout << "  <eigenvalue_sum>  " << setprecision(8)
Francois Gygi committed
944
                 << eigenvalue_sum << " </eigenvalue_sum>" << endl;
945

Francois Gygi committed
946
          tmap["wf_update"].start();
Francois Gygi committed
947
          wf_stepper->update(dwf);
Francois Gygi committed
948
          tmap["wf_update"].stop();
949

950
          // compare delta_eig_sum only after first iteration
951 952
          if ( ite > 0 )
          {
953
            const double delta_ratio = 0.1;
954
            double delta_eig_sum = fabs(eigenvalue_sum - eigenvalue_sum_m);
955 956
            nonscf_converged |= (delta_eig_sum < delta_ratio * delta_ehart);
#if DEBUG
957
            if ( onpe0 )
958 959 960 961
            {
              cout << " BOSampleStepper::step delta_eig_sum: "
                   << delta_eig_sum << endl;
            }
962
#endif
963 964 965
          }
          ite++;
        }
966

Francois Gygi committed
967 968 969
        // subspace diagonalization
        if ( compute_eigvec || s_.ctrl.wf_diag == "EIGVAL" )
        {
Francois Gygi committed
970
          tmap["energy"].start();
971
          ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
972 973
          tmap["energy"].stop();
          tmap["wfdiag"].start();
Francois Gygi committed
974
          s_.wf.diag(dwf,compute_eigvec);
Francois Gygi committed
975
          tmap["wfdiag"].stop();
976
          if ( onpe0 )
977
          {
Francois Gygi committed
978
            cout << "<eigenset>" << endl;
979 980 981 982 983
            // print eigenvalues
            for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
            {
              for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
              {
984 985 986
                const int nst = wf.sd(ispin,ikp)->nst();
                const double eVolt = 2.0 * 13.6058;
                cout <<    "  <eigenvalues spin=\"" << ispin
987 988 989
                     << "\" kpoint=\""
                     << setprecision(8)
                     << wf.sd(ispin,ikp)->kpoint()
990 991 992
                     << "\" weight=\""
                     << setprecision(8)
                     << wf.weight(ikp)
993 994
                     << "\" n=\"" << nst << "\">" << endl;
                for ( int i = 0; i < nst; i++ )
995
                {
996 997 998
                  cout << setw(12) << setprecision(5)
                       << wf.sd(ispin,ikp)->eig(i)*eVolt;
                  if ( i%5 == 4 ) cout << endl;
999
                }
1000 1001
                if ( nst%5 != 0 ) cout << endl;
                cout << "  </eigenvalues>" << endl;
1002
              }
1003
            }
Francois Gygi committed
1004
            cout << "</eigenset>" << endl;
1005
          }
1006
        }
1007

1008
        // update occupation numbers if fractionally occupied states
1009
        // compute weighted sum of eigenvalues
1010 1011 1012
        // default value if no fractional occupation
        double w_eigenvalue_sum = 2.0 * eigenvalue_sum;

1013
        if ( fractional_occ )
Francois Gygi committed
1014
        {
1015
          wf.update_occ(s_.ctrl.fermi_temp);
1016
#if 0
Francois Gygi committed
1017
          if ( onpe0 )
Francois Gygi committed
1018
          {
1019 1020 1021
            cout << "  Wavefunction entropy: " << wf_entropy << endl;
            cout << "  Entropy contribution to free energy: "
                 << - wf_entropy * s_.ctrl.fermi_temp * boltz << endl;
Francois Gygi committed
1022
          }
1023
#endif
1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041
          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() +
1042
                     ef_.esr() - ef_.eself() + ef_.dxc() + ef_.ets();
1043 1044 1045 1046 1047 1048 1049 1050 1051 1052
#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
1053
        }
1054 1055 1056 1057 1058 1059 1060 1061 1062 1063
#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";
        }

1064 1065
        etotal_mm = etotal_m;
        etotal_m = etotal;
1066
        etotal = etotal_int;
1067

1068
        if ( nite_ > 0 && onpe0 )
1069
          cout << "  BOSampleStepper: end scf iteration" << endl;
1070

1071
        // delta_etotal = interval containing etotal, etotal_m and etotal_mm
1072
        double delta_etotal = fabs(etotal - etotal_m);
1073 1074
        delta_etotal = max(delta_etotal,fabs(etotal - etotal_mm));
        delta_etotal = max(delta_etotal,fabs(etotal_m - etotal_mm));
1075 1076 1077
        scf_converged |= (delta_etotal < s_.ctrl.scf_tol);
        itscf++;
      } // while scf
1078

1079 1080
      if ( compute_mlwf || compute_mlwfc )
      {
Francois Gygi committed
1081
        tmap["mlwf"].start();
1082
        for ( int ispin = 0; ispin < nspin; ispin++ )
1083
        {
1084
          mlwft[ispin]->update();
1085
          mlwft[ispin]->compute_transform();
1086
        }
1087 1088

        if ( compute_mlwf )
1089
        {
1090
          for ( int ispin = 0; ispin < nspin; ispin++ )
1091
          {
1092 1093
            SlaterDet& sd = *(wf.sd(ispin,0));
            mlwft[ispin]->apply_transform(sd);
1094 1095
          }
        }
1096 1097 1098

        if ( onpe0 )
        {
1099
          D3vector edipole_sum;
1100
          cout << "<mlwfs>" << endl;
1101
          for ( int ispin = 0; ispin < nspin; ispin++ )
1102
          {
1103
            SlaterDet& sd = *(wf.sd(ispin,0));
1104
            cout << " <mlwfset spin=\"" << ispin
1105
                 << "\" size=\"" << sd.nst() << "\">" << endl;
1106 1107 1108
            double total_spread[6];
            for ( int j = 0; j < 6; j++ )
               total_spread[j] = 0.0;
1109 1110
            for ( int i = 0; i < sd.nst(); i++ )
            {
1111
              D3vector ctr = mlwft[ispin]->center(i);
1112 1113 1114 1115 1116 1117 1118
              double spi[6];
              for (int j=0; j<3; j++)
              {
                spi[j] = mlwft[ispin]->spread2(i,j);
                total_spread[j] += spi[j];
              }

1119 1120 1121 1122 1123 1124
              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
1125 1126 1127 1128
                   << " \" spread=\" "
                   << setw(12) << spi[0]
                   << setw(12) << spi[1]
                   << setw(12) << spi[2] << " \"/>"
1129 1130 1131
                   << endl;
            }

1132 1133 1134 1135
            cout << " <total_spread> ";
            for ( int j = 0; j < 3; j++ )
              cout << setw(10) << total_spread[j];
            cout << " </total_spread>" << endl;
1136 1137
            D3vector edipole = mlwft[ispin]->dipole();
            cout << " <electronic_dipole spin=\"" << ispin << "\"> " << edipole
1138
                 << " </electronic_dipole>" << endl;
1139
            edipole_sum += edipole;
1140
            cout << " </mlwfset>" << endl;
1141
          }
1142 1143 1144 1145 1146 1147 1148
          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;
1149
          cout << "</mlwfs>" << endl;
1150
        }
Francois Gygi committed
1151
        tmap["mlwf"].stop();
1152 1153
      }

1154 1155 1156 1157
      // If GS calculation only, print energy and atomset at end of iterations
      if ( gs_only )
      {
        tmap["charge"].start();
1158
        cd_.update_density();
1159 1160
        tmap["charge"].stop();

Francois Gygi committed
1161
        tmap["update_vhxc"].start();
1162
        ef_.update_vhxc(compute_stress);
Francois Gygi committed
1163
        tmap["update_vhxc"].stop();
1164
        const bool compute_forces = true;
Francois Gygi committed
1165
        tmap["energy"].start();
1166
        ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
Francois Gygi committed
1167
        tmap["energy"].stop();
1168

Francois Gygi committed
1169
        if ( onpe0 )
1170
        {
Francois Gygi committed
1171
          cout << ef_;
1172 1173
          if ( ef_.el_enth() )
            cout << *ef_.el_enth();
1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196
          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
1197 1198 1199 1200 1201 1202 1203
          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
1204
          cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
1205
               << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
1206
          cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
1207
               << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
1208
          cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
1209
               << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
1210
          cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
1211
               << atoms.cell().volume() << " </unit_cell_volume>" << endl;
1212 1213 1214 1215 1216 1217 1218
          if ( compute_stress )
          {
            compute_sigma();
            print_stress();
          }
        }
      }
1219
      wf_stepper->postprocess();
1220
    }
1221 1222 1223 1224
    else
    {
      // wf_stepper == 0, wf_dyn == LOCKED
      // evaluate and print energy
Francois Gygi committed
1225
      tmap["charge"].start();
Francois Gygi committed
1226
      cd_.update_density();
Francois Gygi committed
1227
      tmap["charge"].stop();