BOSampleStepper.C 41.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 "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

        if ( ef_.el_enth() )
          cout << *ef_.el_enth();
406
      }
407

408 409 410 411 412 413 414 415 416 417 418 419
      if ( iter > 0 && ionic_stepper )
      {
        ionic_stepper->compute_v(energy,fion);
      }
      // at this point, positions r0, velocities v0 and forces fion are
      // consistent
      double ekin_ion = 0.0, temp_ion = 0.0;
      if ( ionic_stepper )
      {
        ekin_ion = ionic_stepper->ekin();
        temp_ion = ionic_stepper->temp();
      }
420

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

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

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

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

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

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

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

507
    // Recalculate ground state wavefunctions
508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529
#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
530
    // wavefunction extrapolation
Francois Gygi committed
531
    if ( atoms_move && extrapolate_wf )
532
    {
Francois Gygi committed
533
      for ( int ispin = 0; ispin < nspin; ispin++ )
534 535 536
      {
        for ( int ikp = 0; ikp < s_.wf.nkp(); ikp++ )
        {
537
          if ( ntc_extrapolation )
538
          {
539 540 541 542 543 544 545
            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 )
546
            {
547
              // copy c on cv
548
              for ( int i = 0; i < len; i++ )
549
              {
550 551 552 553 554
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
555
              }
556 557 558
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
559 560 561 562 563
            }
            else if ( iter == 1 )
            {
              s_.wfv->align(s_.wf);
              for ( int i = 0; i < len; i++ )
564
              {
565 566 567 568 569
                const double x = c[i];
                const double xm = cv[i];
                c[i] = 2.0 * x - xm;
                cv[i] = x;
                cmm[i] = xm;
570
              }
571 572 573 574 575 576 577 578 579 580 581 582
              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++ )
583
              {
584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604
                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];
605
              }
606 607 608 609 610 611 612 613 614 615 616 617 618 619 620
            }
            // 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++ )
              {
621 622 623 624 625
                const double x = c[i];
                const double v = cv[i];
                // extrapolation using velocity in cv
                c[i] = x + dt * v;
                cv[i] = x;
626
              }
627 628 629
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
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 684 685 686
            }
            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 )
            {
687
              // copy c to cv
688 689
              for ( int i = 0; i < len; i++ )
              {
690 691 692 693
                const double x = c[i];
                const double v = cv[i];
                c[i] = x + dt * v;
                cv[i] = x;
694
              }
695 696 697 698 699 700
              //tmap["lowdin"].start();
              //s_.wf.sd(ispin,ikp)->lowdin();
              //tmap["lowdin"].stop();
              tmap["gram"].start();
              s_.wf.sd(ispin,ikp)->gram();
              tmap["gram"].stop();
701 702 703
            }
            else
            {
704 705 706 707
              tmap["align"].start();
              s_.wfv->align(s_.wf);
              tmap["align"].stop();

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

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

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

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

743 744
      double ehart, ehart_m;

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

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

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

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

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

792 793 794 795
              const Context * const kpctxt = s_.wf.kpcontext();
              if ( kpctxt->mycol() == 0 )
              {
                // use AndersonMixer on first column only and bcast results
796
                for ( int ispin = 0; ispin < nspin; ispin++ )
797
                {
798 799 800 801 802 803 804
                  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);
805
                }
806 807 808
              }
              else
              {
809
                for ( int ispin = 0; ispin < nspin; ispin++ )
810
                {
811 812 813 814 815
                  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);
816
                }
817
              }
818

819
              for ( int ispin = 0; ispin < nspin; ispin++ )
820
              {
821 822 823 824 825
                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];
826
              }
827 828
            }
            else
829
            {
830
              for ( int ispin = 0; ispin < nspin; ispin++ )
831
              {
832
                for ( int i = 0; i < rhog_in[ispin].size(); i++ )
833
                  rhog_in[ispin][i] += alpha * drhog[ispin][i] * wkerker[i];
834
              }
835
            }
836

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

845
        ef_.update_vhxc();
846

847
        // reset stepper only if multiple non-selfconsistent steps
Francois Gygi committed
848
        if ( nite_ > 1 ) wf_stepper->preprocess();
849

850
        // non-self-consistent loop
851 852 853
        // 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
854
        bool nscf_converged = false;
855 856
        if ( itscf > 0 )
          ehart_m = ehart;
857
        ehart = ef_.ehart();
858 859 860
        double delta_ehart = 0.0;
        if ( itscf > 0 )
          delta_ehart = fabs(ehart - ehart_m);
Francois Gygi committed
861 862
        // if ( onpe0 && nite_ > 1 )
        //   cout << " delta_ehart = " << delta_ehart << endl;
863
        int ite = 0;
864 865
        double etotal_int, etotal_int_m;
        double eigenvalue_sum, eigenvalue_sum_m;
866
        while ( !nscf_converged && ite < nite_ )
Francois Gygi committed
867
        {
868
          double energy = ef_.energy(true,dwf,false,fion,false,sigma_eks);
869
          double enthalpy = energy;
870

871 872 873
          if ( ite > 0 )
            etotal_int_m = etotal_int;

874 875
          etotal_int = energy;

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

883 884 885
          if ( ite > 0 )
            eigenvalue_sum_m = eigenvalue_sum;

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

Francois Gygi committed
891
          wf_stepper->update(dwf);
892

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

          // 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 )
916
            {
917 918
              cout << " BOSampleStepper::step: delta_e_int: "
                   << delta_etotal_int << endl;
919 920 921
              cout << " BOSampleStepper::step: delta_ehart: "
                   << delta_ehart << endl;
            }
922 923 924 925
#else
            double delta_eig_sum = fabs(eigenvalue_sum - eigenvalue_sum_m);
            nscf_converged |= (delta_eig_sum < 0.2 * delta_ehart);
            if ( onpe0 )
926 927 928 929 930 931
            {
              cout << " BOSampleStepper::step delta_eig_sum: "
                   << delta_eig_sum << endl;
              cout << " BOSampleStepper::step: delta_ehart: "
                   << delta_ehart << endl;
            }
932 933 934 935 936
#endif

          }
          ite++;
        }
937 938 939 940

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

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

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

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

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

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

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

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

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

1066 1067
      if ( onpe0 && ef_.el_enth() )
      {
1068
        cout << *ef_.el_enth();
1069 1070
      }

1071 1072 1073 1074
      // If GS calculation only, print energy and atomset at end of iterations
      if ( gs_only )
      {
        tmap["charge"].start();
1075
        cd_.update_density();
1076 1077 1078 1079 1080 1081
        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
1082
        if ( onpe0 )
1083
        {
Francois Gygi committed
1084
          cout << ef_;
1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107
          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
1108 1109 1110 1111 1112 1113 1114
          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
1115
          cout << setprecision(3) << "<unit_cell_alpha>  "
Francois Gygi committed
1116
               << atoms.cell().alpha() << " </unit_cell_alpha>" << endl;
Francois Gygi committed
1117
          cout << setprecision(3) << "<unit_cell_beta>   "
Francois Gygi committed
1118
               << atoms.cell().beta() << " </unit_cell_beta>" << endl;
Francois Gygi committed
1119
          cout << setprecision(3) << "<unit_cell_gamma>  "
Francois Gygi committed
1120
               << atoms.cell().gamma() << " </unit_cell_gamma>" << endl;
Francois Gygi committed
1121
          cout << setprecision(3) << "<unit_cell_volume> "
Francois Gygi committed
1122
               << atoms.cell().volume() << " </unit_cell_volume>" << endl;
1123 1124 1125 1126 1127 1128 1129
          if ( compute_stress )
          {
            compute_sigma();
            print_stress();
          }
        }
      }
1130
      wf_stepper->postprocess();
1131
    }
1132 1133 1134 1135
    else
    {
      // wf_stepper == 0, wf_dyn == LOCKED
      // evaluate and print energy
Francois Gygi committed
1136
      tmap["charge"].start();
Francois Gygi committed
1137
      cd_.update_density();
Francois Gygi committed
1138
      tmap["charge"].stop();
Francois Gygi committed
1139
      ef_.update_vhxc();
Francois Gygi committed
1140
      ef_.energy(true,dwf,false,fion,false,sigma_eks);
Francois Gygi committed
1141
      if ( onpe0 )
1142 1143 1144 1145
      {
        cout << ef_;
      }
    }
1146

Francois Gygi committed
1147 1148 1149
#ifdef USE_APC
    ApcStop(1);
#endif
1150 1151 1152 1153 1154 1155
    // 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
1156
    if ( onpe0 )
1157
    {
1158 1159 1160 1161 1162
      cout << "  <timing name=\"iteration\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
      cout << "</iteration>" << endl;
1163
    }
Francois Gygi committed
1164
    if ( atoms_move )
Francois Gygi committed
1165
      s_.constraints.update_constraints(dt);
1166
  } // for iter
1167

Francois Gygi committed
1168
  if ( atoms_move )
1169
  {
1170 1171
    // compute ionic forces at last position to update velocities
    // consistently with last position
Francois Gygi committed
1172
    tmap["charge"].start();
Francois Gygi committed
1173
    cd_.update_density();
Francois Gygi committed
1174
    tmap["charge"].stop();
1175

Francois Gygi committed
1176
    ef_.update_vhxc();
Francois Gygi committed
1177
    const bool compute_forces = true;
1178
    double energy =
1179
      ef_.energy(false,dwf,compute_forces,fion,compute_stress,sigma_eks);
1180

Francois Gygi committed
1181 1182
    ionic_stepper->compute_v(energy,fion);
    // positions r0 and velocities v0 are consistent
1183
  }
1184

1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251
  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;
  }

1252

1253
  for ( int ispin = 0; ispin < nspin; ispin++ )
1254
  {
1255 1256
   delete mlwft[ispin];
   delete mixer[ispin];
1257
  }
1258

1259
  // delete steppers
1260 1261 1262
  delete wf_stepper;
  delete ionic_stepper;
  delete cell_stepper;
1263

1264
  if ( ntc_extrapolation || asp_extrapolation ) delete wfmm;
1265 1266

  initial_atomic_density = false;
1267
}