BOSampleStepper.C 43.3 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 "SDIonicStepper.h"
Francois Gygi committed
29
#include "SDAIonicStepper.h"
Francois Gygi committed
30
#include "CGIonicStepper.h"
31
#include "MDIonicStepper.h"
Francois Gygi committed
32
#include "BMDIonicStepper.h"
33
#include "SDCellStepper.h"
34
#include "CGCellStepper.h"
35
#include "AndersonMixer.h"
36
#include "MLWFTransform.h"
37
#include "D3tensor.h"
38

Francois Gygi committed
39 40 41 42
#ifdef USE_APC
#include "apc.h"
#endif

43 44 45 46 47
#include <iostream>
#include <iomanip>
using namespace std;

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

Francois Gygi committed
53 54 55 56 57 58 59 60 61 62 63 64
////////////////////////////////////////////////////////////////////////////////
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 )
    {
65 66 67 68 69
      cout << "<timing name=\""
           << setw(15) << (*i).first << "\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
Francois Gygi committed
70 71 72 73
    }
  }
}

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

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

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

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

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

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

  cd_.rhog[0] = rhopst;
  initial_atomic_density = true;
}

153
////////////////////////////////////////////////////////////////////////////////
154
void BOSampleStepper::step(int niter)
155
{
156 157
  const Context& ctxt = s_.ctxt_;
  const bool onpe0 = ctxt.onpe0();
Francois Gygi committed
158

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

161 162 163 164 165
  // 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";
166 167
  const bool compute_mlwf = s_.ctrl.wf_diag == "MLWF";
  const bool compute_mlwfc = s_.ctrl.wf_diag == "MLWFC";
168
  enum ortho_type { GRAM, LOWDIN, ORTHO_ALIGN, RICCATI };
169

170 171
  AtomSet& atoms = s_.atoms;
  Wavefunction& wf = s_.wf;
Francois Gygi committed
172
  const int nspin = wf.nspin();
173

174
  const UnitCell& cell = wf.cell();
175

176
  const double dt = s_.ctrl.dt;
177

178 179
  const string wf_dyn = s_.ctrl.wf_dyn;
  const string atoms_dyn = s_.ctrl.atoms_dyn;
180
  const string cell_dyn = s_.ctrl.cell_dyn;
181

182
  const bool extrapolate_wf = ( atoms_dyn == "MD" );
Francois Gygi committed
183

184
  const bool ntc_extrapolation =
185
    s_.ctrl.debug.find("NTC_EXTRAPOLATION") != string::npos;
186
  const bool asp_extrapolation =
187
    s_.ctrl.debug.find("ASP_EXTRAPOLATION") != string::npos;
188 189

  Wavefunction* wfmm;
Francois Gygi committed
190 191
  if ( extrapolate_wf && ( ntc_extrapolation || asp_extrapolation ) )
    wfmm = new Wavefunction(wf);
192

Francois Gygi committed
193 194
  // Next lines: special value of niter = 0: GS calculation only
  const bool atoms_move = ( niter > 0 && atoms_dyn != "LOCKED" );
195
  const bool compute_stress = ( s_.ctrl.stress == "ON" );
Francois Gygi committed
196 197
  const bool cell_moves = ( niter > 0 && compute_stress &&
                            cell_dyn != "LOCKED" );
198 199
  // GS-only calculation:
  const bool gs_only = !atoms_move && !cell_moves;
200
  const bool use_confinement = ( s_.ctrl.ecuts > 0.0 );
201

202
  Timer tm_iter;
203

204 205
  WavefunctionStepper* wf_stepper = 0;
  if ( wf_dyn == "SD" )
206 207 208
  {
    const double emass = s_.ctrl.emass;
    double dt2bye = (emass == 0.0) ? 0.5 / wf.ecut() : dt*dt/emass;
209

210 211 212 213 214 215 216 217
    // 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);
  }
218
  else if ( wf_dyn == "PSD" )
Francois Gygi committed
219
    wf_stepper = new PSDWavefunctionStepper(wf,s_.ctrl.ecutprec,tmap);
220
  else if ( wf_dyn == "PSDA" )
Francois Gygi committed
221
    wf_stepper = new PSDAWavefunctionStepper(wf,s_.ctrl.ecutprec,tmap);
222
  else if ( wf_dyn == "JD" )
Francois Gygi committed
223
    wf_stepper = new JDWavefunctionStepper(wf,s_.ctrl.ecutprec,ef_,tmap);
224

225
  // wf_stepper == 0 indicates that wf_dyn == LOCKED
226

227 228 229
  IonicStepper* ionic_stepper = 0;
  if ( atoms_dyn == "SD" )
    ionic_stepper = new SDIonicStepper(s_);
Francois Gygi committed
230 231
  else if ( atoms_dyn == "SDA" )
    ionic_stepper = new SDAIonicStepper(s_);
Francois Gygi committed
232 233
  else if ( atoms_dyn == "CG" )
    ionic_stepper = new CGIonicStepper(s_);
234 235
  else if ( atoms_dyn == "MD" )
    ionic_stepper = new MDIonicStepper(s_);
Francois Gygi committed
236 237
  else if ( atoms_dyn == "BMD" )
    ionic_stepper = new BMDIonicStepper(s_);
238

Francois Gygi committed
239 240
  if ( ionic_stepper )
    ionic_stepper->setup_constraints();
241

242 243 244
  CellStepper* cell_stepper = 0;
  if ( cell_dyn == "SD" )
    cell_stepper = new SDCellStepper(s_);
245 246
  else if ( cell_dyn == "CG" )
    cell_stepper = new CGCellStepper(s_);
247

248
  // Allocate wavefunction velocity if not available
Francois Gygi committed
249
  if ( atoms_move && extrapolate_wf )
250 251
  {
    if ( s_.wfv == 0 )
252
    {
253
      s_.wfv = new Wavefunction(wf);
254 255 256 257
      s_.wfv->clear();
    }
  }

258
  vector<MLWFTransform*> mlwft(nspin);
259 260 261 262 263 264 265 266 267 268 269 270 271 272 273

  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;
    }
274 275 276

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

279
  // Charge mixing variables
280 281 282 283 284 285 286 287 288 289 290
  vector<vector<complex<double> > > rhog_in(nspin);
  vector<vector<complex<double> > > drhog(nspin);
  vector<vector<complex<double> > > rhobar(nspin);
  vector<vector<complex<double> > > drhobar(nspin);
  for ( int ispin = 0; ispin < nspin; ispin++ )
  {
    rhog_in[ispin].resize(cd_.rhog[ispin].size());
    drhog[ispin].resize(rhog_in[ispin].size());
    rhobar[ispin].resize(rhog_in[ispin].size());
    drhobar[ispin].resize(rhog_in[ispin].size());
  }
291
  const int anderson_ndim = s_.ctrl.charge_mix_ndim;
292 293 294
  vector<double> wkerker(cd_.rhog[0].size());
  vector<double> wls(cd_.rhog[0].size());

295 296
  vector<AndersonMixer*> mixer(nspin);
  for ( int ispin = 0; ispin < nspin; ispin++ )
297 298
    mixer[ispin] =
      new AndersonMixer(2*rhog_in[ispin].size(),anderson_ndim,&cd_.vcontext());
299 300 301 302 303

  // 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();
304 305 306 307 308 309 310 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;
    cout << " override rc1 value: rc1 = " << rc1 << endl;
    assert(rc1 >= 0.0);
  }

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

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

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

354
    tm_iter.start();
Francois Gygi committed
355 356 357
#ifdef USE_APC
    ApcStart(1);
#endif
358

Francois Gygi committed
359
    if ( onpe0 )
360
      cout << "<iteration count=\"" << iter+1 << "\">\n";
361

Francois Gygi committed
362
    // compute energy and ionic forces using existing wavefunction
Francois Gygi committed
363

364 365 366
    if ( !gs_only )
    {
      tmap["charge"].start();
367
      cd_.update_density();
368
      tmap["charge"].stop();
369

370 371 372 373
      ef_.update_vhxc();
      const bool compute_forces = true;
      double energy =
        ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
374
      double enthalpy = energy;
375

376
      if ( onpe0 )
377
      {
378 379 380 381 382 383 384 385 386 387 388 389
        cout.setf(ios::fixed,ios::floatfield);
        cout.setf(ios::right,ios::adjustfield);
        cout << "  <ekin>   " << setprecision(8)
             << setw(15) << ef_.ekin() << " </ekin>\n";
        if ( use_confinement )
          cout << "  <econf>  " << setw(15) << ef_.econf() << " </econf>\n";
        cout << "  <eps>    " << setw(15) << ef_.eps() << " </eps>\n"
             << "  <enl>    " << setw(15) << ef_.enl() << " </enl>\n"
             << "  <ecoul>  " << setw(15) << ef_.ecoul() << " </ecoul>\n"
             << "  <exc>    " << setw(15) << ef_.exc() << " </exc>\n"
             << "  <esr>    " << setw(15) << ef_.esr() << " </esr>\n"
             << "  <eself>  " << setw(15) << ef_.eself() << " </eself>\n"
390 391 392 393
             << "  <ets>    " << setw(15) << ef_.ets() << " </ets>\n";
        if ( s_.extforces.size() > 0 )
          cout << "  <eexf>     " << setw(15) << ef_.eexf() << " </eexf>\n";
        cout << "  <etotal> " << setw(15) << ef_.etotal() << " </etotal>\n";
394 395 396
        if ( compute_stress )
        {
          const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
397
          enthalpy = ef_.etotal() + pext * cell.volume();
398 399 400 401 402
          cout << "  <pv>     " << setw(15) << pext * cell.volume()
               << " </pv>" << endl;
          cout << "  <enthalpy> " << setw(15) << enthalpy << " </enthalpy>\n"
             << flush;
        }
403
      }
404

405 406 407 408 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
      double ekin_ion = 0.0, temp_ion = 0.0;
      if ( ionic_stepper )
      {
        ekin_ion = ionic_stepper->ekin();
        temp_ion = ionic_stepper->temp();
      }
417

418 419
      // print positions, velocities and forces at time t0
      if ( onpe0 )
420
      {
421 422 423
        cout << "<atomset>" << endl;
        cout << atoms.cell();
        for ( int is = 0; is < atoms.atom_list.size(); is++ )
424
        {
425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441
          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
442
        }
443
        cout << "</atomset>" << endl;
Francois Gygi committed
444 445 446 447 448 449 450
        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
451
        cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
452
             << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
453
        cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
454
             << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
455
        cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
456
             << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
457
        cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
458
             << atoms.cell().volume() << " </unit_cell_volume>" << endl;
459 460 461 462 463 464

        // 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
465
        cout << setprecision(8);
466
        cout << "  <econst> " << energy+ekin_ion+ekin_stepper << " </econst>\n";
467 468
        cout << "  <ekin_ion> " << ekin_ion << " </ekin_ion>\n";
        cout << "  <temp_ion> " << temp_ion << " </temp_ion>\n";
Francois Gygi committed
469
      }
470

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

486
      if ( compute_stress )
487
      {
488 489
        compute_sigma();
        print_stress();
490

491 492
        if ( cell_moves )
        {
493
          cell_stepper->compute_new_cell(enthalpy,sigma,fion);
494 495 496

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

498 499 500
          ef_.cell_moved();
          ef_.atoms_moved(); // modifications of the cell also move ions
        }
501
      }
502
    } // if !gs_only
503

504
    // Recalculate ground state wavefunctions
505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
#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
527
    // wavefunction extrapolation
Francois Gygi committed
528
    if ( atoms_move && extrapolate_wf )
529
    {
Francois Gygi committed
530
      for ( int ispin = 0; ispin < nspin; ispin++ )
531 532 533
      {
        for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
        {
534
          if ( ntc_extrapolation )
535
          {
536 537 538 539 540 541 542
            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 )
543
            {
544
              // copy c on cv
545
              for ( int i = 0; i < len; i++ )
546
              {
547 548 549 550 551
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
552
              }
553 554 555
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
556 557 558 559 560
            }
            else if ( iter == 1 )
            {
              s_.wfv->align(s_.wf);
              for ( int i = 0; i < len; i++ )
561
              {
562 563 564 565 566
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
                cmm[i] = xm;
567
              }
568 569 570 571 572 573 574 575 576 577 578 579
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
            }
            else
            {
              // align wf with wfmm before extrapolation
              // s_.wf.align(*wfmm);
              wfmm->align(s_.wf);

              // extrapolate
              for ( int i = 0; i < len; i++ )
580
              {
581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601
                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];
602
              }
603 604 605 606 607 608 609 610 611 612 613 614 615 616 617
            }
            // 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++ )
              {
618 619 620 621 622
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
623
              }
624 625 626
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
            }
            else if ( iter == 1 )
            {
              //s_.wfv->align(s_.wf);
              for ( int i = 0; i < len; i++ )
              {
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
                cmm[i] = xm;
              }
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
            }
            else
            {
              // align wf with wfmm before extrapolation
              // s_.wf.align(*wfmm);
              // wfmm->align(s_.wf);

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

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

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

705 706 707 708 709 710 711 712
              // 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;
              }
713 714 715
              //tmap["ortho_align"].start();
              //s_.wf.sd(ispin,ikp)->ortho_align(*s_.wfv->sd(ispin,ikp));
              //tmap["ortho_align"].stop();
716 717 718 719

              //tmap["riccati"].start();
              //s_.wf.sd(ispin,ikp)->riccati(*s_.wfv->sd(ispin,ikp));
              //tmap["riccati"].stop();
720 721 722 723

              tmap["lowdin"].start();
              s_.wf.sd(ispin,ikp)->lowdin();
              tmap["lowdin"].stop();
724
            }
725 726 727
          }
        }
      }
Francois Gygi committed
728
    } // atoms_move && extrapolate_wf
729

Francois Gygi committed
730
    // do nitscf self-consistent iterations, each with nite electronic steps
731
    if ( wf_stepper != 0 )
732
    {
733
      wf_stepper->preprocess();
Francois Gygi committed
734
      if ( anderson_charge_mixing )
735
      {
736 737
        for ( int ispin = 0; ispin < nspin; ispin++ )
          mixer[ispin]->restart();
738
      }
Francois Gygi committed
739

740 741
      double ehart, ehart_m;

Francois Gygi committed
742
      for ( int itscf = 0; itscf < nitscf_; itscf++ )
743
      {
Francois Gygi committed
744
        if ( nite_ > 1 && onpe0 )
745
          cout << "  BOSampleStepper: start scf iteration" << endl;
746

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

Francois Gygi committed
755
        // charge mixing
756
        if ( nite_ > 1 )
757
        {
758
          if ( itscf == 0 )
Francois Gygi committed
759
          {
760
            // at first scf iteration, nothing to mix
761
            // memorize rhog in rhog_in for next iteration
762

763 764
            for ( int ispin = 0; ispin < nspin; ispin++ )
              rhog_in[ispin] = cd_.rhog[ispin];
765 766 767
          }
          else
          {
768 769
            // itscf > 0
            // compute unscreened correction drhog
770
            for ( int ispin = 0; ispin < nspin; ispin++ )
771
            {
772
              for ( int i = 0; i < rhog_in[ispin].size(); i++ )
773
              {
774
                drhog[ispin][i] = (cd_.rhog[ispin][i] - rhog_in[ispin][i]);
775 776
              }
            }
777 778 779 780

            const double alpha = s_.ctrl.charge_mix_coeff;
            // Anderson acceleration
            if ( anderson_charge_mixing )
781
            {
782
              // row weighting of LS calculation
783
              for ( int ispin = 0; ispin < nspin; ispin++ )
784
              {
785 786
                for ( int i = 0; i < drhog[ispin].size(); i++ )
                  drhog[ispin][i] /= wls[i];
787
              }
788

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

816
              for ( int ispin = 0; ispin < nspin; ispin++ )
817
              {
818 819 820 821 822
                for ( int i = 0; i < drhog[ispin].size(); i++ )
                  drhobar[ispin][i] *= wls[i];
                for ( int i = 0; i < rhog_in[ispin].size(); i++ )
                  rhog_in[ispin][i] = rhobar[ispin][i] +
                    alpha * drhobar[ispin][i] * wkerker[i];
823
              }
824 825
            }
            else
826
            {
827
              for ( int ispin = 0; ispin < nspin; ispin++ )
828
              {
829
                for ( int i = 0; i < rhog_in[ispin].size(); i++ )
830
                  rhog_in[ispin][i] += alpha * drhog[ispin][i] * wkerker[i];
831
              }
832
            }
833

834
            for ( int ispin = 0; ispin < nspin; ispin++ )
835
            {
836
              cd_.rhog[ispin] = rhog_in[ispin];
837
            }
838
            cd_.update_rhor();
839
          }
840
        } // if nite_ > 1
Francois Gygi committed
841

842
        ef_.update_vhxc();
843

844
        // reset stepper only if multiple non-selfconsistent steps
Francois Gygi committed
845
        if ( nite_ > 1 ) wf_stepper->preprocess();
846

847
        // non-self-consistent loop
848 849 850
        // repeat until the change in etotal_int or in the
        // eigenvalue sum is smaller than a fraction of the change in
        // Hartree energy in the last scf iteration
851
        bool nscf_converged = false;
852 853
        if ( itscf > 0 )
          ehart_m = ehart;
854
        ehart = ef_.ehart();
855 856 857
        double delta_ehart = 0.0;
        if ( itscf > 0 )
          delta_ehart = fabs(ehart - ehart_m);
Francois Gygi committed
858 859
        // if ( onpe0 && nite_ > 1 )
        //   cout << " delta_ehart = " << delta_ehart << endl;
860
        int ite = 0;
861 862
        double etotal_int, etotal_int_m;
        double eigenvalue_sum, eigenvalue_sum_m;
863
        while ( !nscf_converged && ite < nite_ )
Francois Gygi committed
864
        {
865
          double energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
866
          double enthalpy = energy;
867

868 869 870
          if ( ite > 0 )
            etotal_int_m = etotal_int;

871 872
          etotal_int = energy;

Francois Gygi committed
873 874 875 876
          // 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
877 878
          // Note: since the hamiltonian is hermitian and dwf=H*wf
          // the dot product in the following line is real
879

880 881 882
          if ( ite > 0 )
            eigenvalue_sum_m = eigenvalue_sum;

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

Francois Gygi committed
888
          wf_stepper->update(dwf);
889

Francois Gygi committed
890
          if ( onpe0 )
Francois Gygi committed
891 892 893
          {
            cout.setf(ios::fixed,ios::floatfield);
            cout.setf(ios::right,ios::adjustfield);
Francois Gygi committed
894
            cout << "  <etotal_int> " << setprecision(8) << setw(15)
895
                 << energy << " </etotal_int>\n";
Francois Gygi committed
896 897 898
            if ( compute_stress )
            {
              const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
899
              enthalpy = energy + pext * cell.volume();
900
              cout << "  <enthalpy_int> " << setw(15)
Francois Gygi committed
901 902 903
                   << enthalpy << " </enthalpy_int>\n"
                   << flush;
            }
904
          }
905 906 907 908 909 910 911 912

          // compare delta_etotal_int only after first iteration
          if ( ite > 0 )
          {
#if 1
            double delta_etotal_int = fabs(etotal_int - etotal_int_m);
            nscf_converged |= (delta_etotal_int < 0.2 * delta_ehart);
            if ( onpe0 )
913
            {
914 915
              cout << " BOSampleStepper::step: delta_e_int: "
                   << delta_etotal_int << endl;
916 917 918
              cout << " BOSampleStepper::step: delta_ehart: "
                   << delta_ehart << endl;
            }
919 920 921 922
#else
            double delta_eig_sum = fabs(eigenvalue_sum - eigenvalue_sum_m);
            nscf_converged |= (delta_eig_sum < 0.2 * delta_ehart);
            if ( onpe0 )
923 924 925 926 927 928
            {
              cout << " BOSampleStepper::step delta_eig_sum: "
                   << delta_eig_sum << endl;
              cout << " BOSampleStepper::step: delta_ehart: "
                   << delta_ehart << endl;
            }
929 930 931 932 933
#endif

          }
          ite++;
        }
934 935 936 937

        // if ( onpe0 && nite_ > 1 && ite >= nite_ )
        //   cout << " BOSampleStepper::step: nscf loop not converged after "
        //        << nite_ << " iterations" << endl;
938

Francois Gygi committed
939 940 941
        // subspace diagonalization
        if ( compute_eigvec || s_.ctrl.wf_diag == "EIGVAL" )
        {
942
          ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
943
          s_.wf.diag(dwf,compute_eigvec);
944
          if ( onpe0 )
945
          {
Francois Gygi committed
946
            cout << "<eigenset>" << endl;
947 948 949 950 951
            // print eigenvalues
            for ( int ispin = 0; ispin < wf.nspin(); ispin++ )
            {
              for ( int ikp = 0; ikp < wf.nkp(); ikp++ )
              {
952 953 954
                const int nst = wf.sd(ispin,ikp)->nst();
                const double eVolt = 2.0 * 13.6058;
                cout <<    "  <eigenvalues spin=\"" << ispin
955 956 957
                     << "\" kpoint=\""
                     << setprecision(8)
                     << wf.sd(ispin,ikp)->kpoint()
958 959
                     << "\" n=\"" << nst << "\">" << endl;
                for ( int i = 0; i < nst; i++ )
960
                {
961 962 963
                  cout << setw(12) << setprecision(5)
                       << wf.sd(ispin,ikp)->eig(i)*eVolt;
                  if ( i%5 == 4 ) cout << endl;
964
                }
965 966
                if ( nst%5 != 0 ) cout << endl;
                cout << "  </eigenvalues>" << endl;
967
              }
968
            }
Francois Gygi committed
969
            cout << "</eigenset>" << endl;
970
          }
971
        }
972

973 974
        // update occupation numbers if fractionally occupied states
        if ( fractional_occ )
Francois Gygi committed
975
        {
976 977
          wf.update_occ(s_.ctrl.fermi_temp);
          const double wf_entropy = wf.entropy();
Francois Gygi committed
978
          if ( onpe0 )
Francois Gygi committed
979
          {
980
            cout << "  Wavefunction entropy: " << wf_entropy << endl;
Francois Gygi committed
981
            const double boltz = 1.0 / ( 11605.0 * 2.0 * 13.6058 );
982 983
            cout << "  Entropy contribution to free energy: "
                 << - wf_entropy * s_.ctrl.fermi_temp * boltz << endl;
Francois Gygi committed
984 985
          }
        }
986

Francois Gygi committed
987
        if ( nite_ > 1 && onpe0 )
988
          cout << "  BOSampleStepper: end scf iteration" << endl;
Francois Gygi committed
989
      } // for itscf
990

991 992
      if ( compute_mlwf || compute_mlwfc )
      {
993
        for ( int ispin = 0; ispin < nspin; ispin++ )
994
        {
995
          mlwft[ispin]->update();
996
          mlwft[ispin]->compute_transform();
997
        }
998 999

        if ( compute_mlwf )
1000
        {
1001
          for ( int ispin = 0; ispin < nspin; ispin++ )
1002
          {
1003 1004
            SlaterDet& sd = *(wf.sd(ispin,0));
            mlwft[ispin]->apply_transform(sd);
1005 1006
          }
        }
1007 1008 1009

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

1030 1031 1032 1033 1034 1035
              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
1036 1037 1038 1039
                   << " \" spread=\" "
                   << setw(12) << spi[0]
                   << setw(12) << spi[1]
                   << setw(12) << spi[2] << " \"/>"
1040 1041 1042 1043
                   << endl;
            }

            cout << " </mlwf_set>" << endl;
1044 1045 1046 1047
            cout << " <total_spread> ";
            for ( int j = 0; j < 3; j++ )
              cout << setw(10) << total_spread[j];
            cout << " </total_spread>" << endl;
1048 1049
            D3vector edipole = mlwft[ispin]->dipole();
            cout << " <electronic_dipole spin=\"" << ispin << "\"> " << edipole
1050
                 << " </electronic_dipole>" << endl;
1051
            edipole_sum += edipole;
1052
          }
1053 1054 1055 1056 1057 1058 1059
          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;
1060 1061 1062
        }
      }

1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108
      if ( onpe0 && ef_.el_enth() )
      {
        vector<D3vector>& mlwfc = ef_.el_enth()->mlwfc();
        vector<double>& mlwfs = ef_.el_enth()->mlwfs();
        vector<D3vector>& cor_real = ef_.el_enth()->cor_real();
        vector<D3tensor>& quad = ef_.el_enth()->quad();

        int nst = wf.sd(0,0)->nst();

        cout << " <mlwf_set size=\"" << nst << "\">" << endl;
        for ( int i = 0; i < nst; i++ )
        {
          cout.setf(ios::fixed, ios::floatfield);
          cout.setf(ios::right, ios::adjustfield);
          cout << "   <mlwf center=\"" << setprecision(10)
               << setw(16) << mlwfc[i].x
               << setw(16) << mlwfc[i].y
               << setw(16) << mlwfc[i].z
               << " \" spread=\" " << mlwfs[i] << " \"/>" << endl
               << "    <correction_real center=\"" << setprecision(10)
               << setw(16) << cor_real[i].x
               << setw(16) << cor_real[i].y
               << setw(16) << cor_real[i].z
               << " \"/>" << endl
               << "    <mlwf_ref center=\"" << setprecision(10)
               << setw(16) << mlwfc[i].x + cor_real[i].x
               << setw(16) << mlwfc[i].y + cor_real[i].y
               << setw(16) << mlwfc[i].z + cor_real[i].z;

          cout << " \" spread=\" "
               << sqrt ( quad[i].trace() );

          cout << " \"/>" << endl;

          cout << "    <quad>"
               << setw(16) << quad[i][0]
               << setw(16) << quad[i][4]
               << setw(16) << quad[i][8]
               << setw(16) << quad[i][1]
               << setw(16) << quad[i][2]
               << setw(16) << quad[i][5]
               << " </quad>" << endl;
        }
        cout << " </mlwf_set>" << endl;
      }

1109 1110 1111 1112
      // If GS calculation only, print energy and atomset at end of iterations
      if ( gs_only )
      {
        tmap["charge"].start();
1113
        cd_.update_density();
1114 1115 1116 1117 1118 1119
        tmap["charge"].stop();

        ef_.update_vhxc();
        const bool compute_forces = true;
        ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);

Francois Gygi committed
1120
        if ( onpe0 )
1121
        {
Francois Gygi committed
1122
          cout << ef_;
1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145
          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
1146 1147 1148 1149 1150 1151 1152
          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
1153
          cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
1154
               << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
1155
          cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
1156
               << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
1157
          cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
1158
               << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
1159
          cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
1160
               << atoms.cell().volume() << " </unit_cell_volume>" << endl;
1161 1162 1163 1164 1165 1166 1167
          if ( compute_stress )
          {
            compute_sigma();
            print_stress();
          }
        }
      }
1168
      wf_stepper->postprocess();
1169
    }
1170 1171 1172 1173
    else
    {
      // wf_stepper == 0, wf_dyn == LOCKED
      // evaluate and print energy
Francois Gygi committed
1174
      tmap["charge"].start();
Francois Gygi committed
1175
      cd_.update_density();
Francois Gygi committed
1176
      tmap["charge"].stop();
Francois Gygi committed
1177
      ef_.update_vhxc();
Francois Gygi committed
1178
      ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
1179
      if ( onpe0 )
1180 1181 1182 1183
      {
        cout << ef_;
      }
    }
1184

Francois Gygi committed
1185 1186 1187
#ifdef USE_APC
    ApcStop(1);
#endif
1188 1189 1190 1191 1192 1193
    // print iteration time
    double time = tm_iter.real();
    double tmin = time;
    double tmax = time;
    s_.ctxt_.dmin(1,1,&tmin,1);
    s_.ctxt_.dmax(1,1,&tmax,1);
Francois Gygi committed
1194
    if ( onpe0 )
1195
    {
1196 1197 1198 1199 1200
      cout << "  <timing name=\"iteration\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
      cout << "</iteration>" << endl;
1201
    }
Francois Gygi committed
1202
    if ( atoms_move )
Francois Gygi committed
1203
      s_.constraints.update_constraints(dt);
1204
  } // for iter
1205

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

Francois Gygi committed
1214
    ef_.update_vhxc();
Francois Gygi committed
1215
    const bool compute_forces = true;
1216
    double energy =
1217
      ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
1218

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

1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289
  if ( atoms_move && extrapolate_wf )
  {
    // compute wavefunction velocity after last iteration
    // s_.wfv contains the previous wavefunction

    s_.wfv->align(s_.wf);

    for ( int ispin = 0; ispin < s_.wf.nspin(); ispin++ )
    {
      for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
      {
        double* c = (double*) s_.wf.sd(ispin,ikp)->c().cvalptr();
        double* cm = (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;
        const double dt_inv = 1.0 / dt;
        if ( ntc_extrapolation )
        {
          double* cmm = (double*) wfmm->sd(ispin,ikp)->c().cvalptr();
          for ( int i = 0; i < len; i++ )
          {
            const double x = c[i];
            const double xmm = cmm[i];
            cm[i] = dt_inv * ( x - xmm );
          }
          tmap["gram"].start();
          s_.wf.sd(ispin,ikp)->gram();
          tmap["gram"].stop();
        }
        else // normal extrapolation or asp_extrapolation
        {
          for ( int i = 0; i < len; i++ )
          {
            const double x = c[i];
            const double xm = cm[i];
            cm[i] = dt_inv * ( x - xm );
          }
          tmap["gram"].start();
          s_.wf.sd(ispin,ikp)->gram();
          tmap["gram"].stop();
        }
      }
    }

    // compute ionic forces at last position to update velocities
    // consistently with last position
    tmap["charge"].start();
    cd_.update_density();
    tmap["charge"].stop();

    ef_.update_vhxc();
    const bool compute_forces = true;
    double energy =
      ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);

    ionic_stepper->compute_v(energy,fion);
    // positions r0 and velocities v0 are consistent
  }
  else
  {
    // delete wavefunction velocities
    if ( s_.wfv != 0 )
      delete s_.wfv;
    s_.wfv = 0;
  }

1290

1291
  for ( int ispin = 0; ispin < nspin; ispin++ )
1292
  {
1293 1294
   delete mlwft[ispin];
   delete mixer[ispin];
1295
  }
1296

1297
  // delete steppers
1298 1299 1300
  delete wf_stepper;
  delete ionic_stepper;
  delete cell_stepper;
1301

1302
  if ( ntc_extrapolation || asp_extrapolation ) delete wfmm;
1303 1304

  initial_atomic_density = false;