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

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

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

////////////////////////////////////////////////////////////////////////////////
46 47
BOSampleStepper::BOSampleStepper(Sample& s, int nitscf, int nite) :
  SampleStepper(s), cd_(s.wf), ef_(s,cd_),
48
  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
  AtomSet& atoms = s_.atoms;
  Wavefunction& wf = s_.wf;
Francois Gygi committed
183
  const int nspin = wf.nspin();
184

185
  const double dt = s_.ctrl.dt;
186

187 188
  const string wf_dyn = s_.ctrl.wf_dyn;
  const string atoms_dyn = s_.ctrl.atoms_dyn;
189
  const string cell_dyn = s_.ctrl.cell_dyn;
190

191
  const bool extrapolate_wf = ( atoms_dyn == "MD" );
Francois Gygi committed
192

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

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

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

210 211 212
  const double force_tol = s_.ctrl.force_tol;
  const double stress_tol = s_.ctrl.stress_tol;

213
  Timer tm_iter;
214

215 216
  Preconditioner prec(wf,ef_,s_.ctrl.ecutprec);

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

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

238
  // wf_stepper == 0 indicates that wf_dyn == LOCKED
239

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

Francois Gygi committed
252 253
  if ( ionic_stepper )
    ionic_stepper->setup_constraints();
254

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

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

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

  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;
    }
287 288 289

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

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

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

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

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

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

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

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

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

364
    tm_iter.start();
365

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

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

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

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

389 390 391 392 393 394 395 396 397 398 399 400 401 402
      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 )
      {
        for ( int i = 0; i < sigma.size(); i++ )
          maxstress = max(maxstress, fabs(sigma[i]));
      }

403
      if ( onpe0 )
404
      {
405
        cout << cd_;
406
        cout << ef_;
407 408
        if ( ef_.el_enth() )
          cout << *ef_.el_enth();
409
      }
410

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

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

441 442 443 444 445 446
      double ekin_ion = 0.0, temp_ion = 0.0;
      if ( ionic_stepper )
      {
        ekin_ion = ionic_stepper->ekin();
        temp_ion = ionic_stepper->temp();
      }
447

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

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

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

517
      if ( compute_stress )
518
      {
519 520
        compute_sigma();
        print_stress();
521

522 523
        if ( cell_moves )
        {
524
          cell_stepper->compute_new_cell(enthalpy,sigma,fion);
525 526 527

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

529 530 531
          ef_.cell_moved();
          ef_.atoms_moved(); // modifications of the cell also move ions
        }
532
      }
533
    } // if !gs_only
534

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

              // extrapolate
              for ( int i = 0; i < len; i++ )
610
              {
611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631
                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];
632
              }
633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
            }
            // 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++ )
              {
648 649 650 651 652
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
653
              }
654 655 656
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
657 658 659
            }
            else if ( iter == 1 )
            {
660
              s_.wfv->sd(ispin,ikp)->align(*s_.wf.sd(ispin,ikp));
661 662 663 664 665 666 667 668 669 670 671 672 673 674 675
              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
676
              s_.wf.sd(ispin,ikp)->align(*wfmm->sd(ispin,ikp));
677 678 679 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

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

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

              //tmap["riccati"].start();
              //s_.wf.sd(ispin,ikp)->riccati(*s_.wfv->sd(ispin,ikp));
              //tmap["riccati"].stop();
749 750 751 752

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

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

766
      double ehart, ehart_m;
767
      double delta_ehart;
768 769
      bool scf_converged = false;
      int itscf = 0;
770
      double etotal = 0.0, etotal_m = 0.0, etotal_mm = 0.0;
771

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

777
        // update charge density
Francois Gygi committed
778
        tmap["charge"].start();
779 780 781 782
        // 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 ) )
783
          cd_.update_density();
Francois Gygi committed
784
        tmap["charge"].stop();
785

786 787 788
        if ( onpe0 )
          cout << cd_;

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

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

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

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

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

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

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

880
        // reset stepper only if multiple non-selfconsistent steps
881
        if ( nite_ > 0 ) wf_stepper->preprocess();
882

883
        // non-self-consistent loop
884
        // repeat until the change in eigenvalue_sum is smaller than a
885 886
        // fraction delta_ratio of the change in Hartree energy delta_ehart
        // in the last scf iteration
887
        bool nonscf_converged = false;
888 889 890 891 892 893 894
        if ( itscf == 0 )
        {
          ehart = ef_.ehart();
          delta_ehart = 0.0;
        }
        else if ( itscf == 1 )
        {
895
          ehart_m = ehart;
896
          ehart = ef_.ehart();
897
          delta_ehart = fabs(ehart - ehart_m);
898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914
        }
        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

915
        int ite = 0;
Francois Gygi committed
916
        double etotal_int;
917 918 919 920 921 922

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

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

935 936 937
          if ( ite > 0 )
            eigenvalue_sum_m = eigenvalue_sum;

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

Francois Gygi committed
943
          tmap["wf_update"].start();
Francois Gygi committed
944
          wf_stepper->update(dwf);
Francois Gygi committed
945
          tmap["wf_update"].stop();
946

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

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

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

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

1061 1062
        etotal_mm = etotal_m;
        etotal_m = etotal;
1063
        etotal = etotal_int;
1064

1065
        if ( nite_ > 0 && onpe0 )
1066
          cout << "  BOSampleStepper: end scf iteration" << endl;
1067

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

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

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

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

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

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

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

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

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