CPSampleStepper.C 8.51 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
3
// CPSampleStepper.C
4 5
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: CPSampleStepper.C,v 1.19 2008-04-15 01:36:15 fgygi Exp $
7 8 9 10 11

#include "CPSampleStepper.h"
#include "SlaterDet.h"
#include "MDWavefunctionStepper.h"
#include "MDIonicStepper.h"
12
#include "SDCellStepper.h"
13
#include "Basis.h"
Francois Gygi committed
14
#include "Species.h"
15 16 17 18 19 20

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

////////////////////////////////////////////////////////////////////////////////
21
CPSampleStepper::CPSampleStepper(Sample& s) :
Francois Gygi committed
22
  SampleStepper(s), cd_(s.wf), ef_(s,cd_), dwf(s.wf), wfv(s.wfv)
23
{
24 25
  const double emass = s.ctrl.emass;
  const double dt = s.ctrl.dt;
26 27 28 29 30 31 32 33
  double dt2bye = (emass == 0.0) ? 0.5 / s.wf.ecut() : dt*dt/emass;

  // divide dt2bye by facs coefficient if stress == ON
  const double facs = 2.0;
  if ( s.ctrl.stress == "ON" )
  {
    dt2bye /= facs;
  }
34 35 36 37 38
  if ( s.wfv == 0 )
  {
    s.wfv = new Wavefunction(s.wf);
    s.wfv->clear();
  }
39
  mdwf_stepper = new MDWavefunctionStepper(s.wf,s.wfv,dt,dt2bye,tmap);
40 41 42 43 44 45 46
  assert(mdwf_stepper!=0);
  mdionic_stepper = 0;
  if ( s.ctrl.atoms_dyn != "LOCKED" )
    mdionic_stepper = new MDIonicStepper(s);
}

////////////////////////////////////////////////////////////////////////////////
47 48 49 50
CPSampleStepper::~CPSampleStepper(void)
{
  delete mdwf_stepper;
  if ( mdionic_stepper != 0 ) delete mdionic_stepper;
Francois Gygi committed
51 52 53 54 55 56 57 58 59
  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 )
    {
60 61 62 63 64
      cout << "<timing name=\""
           << setw(15) << (*i).first << "\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
Francois Gygi committed
65 66
    }
  }
67 68 69 70
}

////////////////////////////////////////////////////////////////////////////////
void CPSampleStepper::step(int niter)
71
{
Francois Gygi committed
72
  const bool onpe0 = s_.ctxt_.onpe0();
73 74 75 76 77 78 79 80 81 82 83 84
  // CP dynamics is allowed only for all doubly occupied states
  // check if states are all doubly occupied
  const bool wf_double_occ = (s_.wf.nel() == 2 * s_.wf.nst());
  if ( !wf_double_occ )
  {
    if ( s_.ctxt_.onpe0() )
    {
      cout << " CPSampleStepper::step:"
              " not all states doubly occupied: cannot run" << endl;
    }
    return;
  }
85

86 87
  AtomSet& atoms = s_.atoms;
  Wavefunction& wf = s_.wf;
88

89
  const double dt = s_.ctrl.dt;
90
  double ekin_ion=0.0,ekin_e, temp_ion=0.0, eta;
91

92 93 94
  const string wf_dyn = s_.ctrl.wf_dyn;
  assert(wf_dyn=="MD");
  const string atoms_dyn = s_.ctrl.atoms_dyn;
95
  const string cell_dyn = s_.ctrl.cell_dyn;
96

97 98
  const bool compute_hpsi = true;
  const bool compute_forces = ( atoms_dyn != "LOCKED" );
99
  const bool compute_stress = ( s_.ctrl.stress == "ON" );
100
  const bool use_confinement = ( s_.ctrl.ecuts > 0.0 );
101

102 103 104
  CellStepper* cell_stepper = 0;
  if ( cell_dyn == "SD" )
    cell_stepper = new SDCellStepper(s_);
105

106 107 108 109 110
  if ( s_.wfv == 0 )
  {
    s_.wfv = new Wavefunction(wf);
    s_.wfv->clear();
  }
111

Francois Gygi committed
112 113
  if ( mdionic_stepper )
    mdionic_stepper->setup_constraints();
114

115
  Timer tm_iter;
116

Francois Gygi committed
117
  tmap["charge"].start();
Francois Gygi committed
118
  cd_.update_density();
Francois Gygi committed
119
  tmap["charge"].stop();
120

Francois Gygi committed
121
  ef_.update_vhxc();
122
  double energy =
123
    ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
124

125 126
  mdwf_stepper->compute_wfm(dwf);

127 128 129 130 131
  for ( int iter = 0; iter < niter; iter++ )
  {
    tm_iter.reset();
    tm_iter.start();
    if ( s_.ctxt_.mype() == 0 )
132
      cout << "<iteration count=\"" << iter+1 << "\">\n";
133

Francois Gygi committed
134 135
    if ( mdionic_stepper )
      atoms.sync();
136

137
    mdwf_stepper->update(dwf);
138

139
    ekin_e = mdwf_stepper->ekin();
140

Francois Gygi committed
141
    if ( onpe0 )
142 143 144
    {
      cout.setf(ios::fixed,ios::floatfield);
      cout.setf(ios::right,ios::adjustfield);
145
      cout << "  <ekin>     " << setprecision(8)
146
           << setw(15) << ef_.ekin() << " </ekin>\n";
147
      if ( use_confinement )
148
      {
149 150
        cout << "  <econf>    " << setw(15) << ef_.econf()
             << " </econf>\n";
151
      }
152 153 154 155 156 157 158
      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"
           << "  <etotal>   " << setw(15) << ef_.etotal() << " </etotal>\n"
159
           << flush;
160 161 162 163 164 165 166 167 168
      if ( compute_stress )
      {
        const double pext = (sigma_ext[0]+sigma_ext[1]+sigma_ext[2])/3.0;
        const double enthalpy = ef_.etotal() + pext * s_.wf.cell().volume();
        cout << "  <pv>     " << setw(15) << pext * s_.wf.cell().volume()
             << " </pv>" << endl;
        cout << "  <enthalpy> " << setw(15) << enthalpy << " </enthalpy>\n"
           << flush;
      }
169
    }
170

171 172
    if ( compute_forces )
    {
173
      if ( iter > 0 )
174
      {
Francois Gygi committed
175
        mdionic_stepper->compute_v(energy,fion);
176
      }
Francois Gygi committed
177
      mdionic_stepper->compute_r(energy,fion);
178
      ekin_ion = mdionic_stepper->ekin();
179

Francois Gygi committed
180
      if ( onpe0 )
181
      {
182
        cout << "<atomset>" << endl;
Francois Gygi committed
183
        cout << atoms.cell();
184 185 186 187 188 189 190 191 192
        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"
193
                 << "    <position> "
194
                 << mdionic_stepper->r0(is,i) << " "
195
                 << mdionic_stepper->r0(is,i+1) << " "
196
                 << mdionic_stepper->r0(is,i+2) << " </position>\n"
197
                 << "    <velocity> "
198
                 << mdionic_stepper->v0(is,i) << " "
199
                 << mdionic_stepper->v0(is,i+1) << " "
200
                 << mdionic_stepper->v0(is,i+2) << " </velocity>\n"
201
                 << "    <force> "
202
                 << fion[is][i] << " "
203
                 << fion[is][i+1] << " "
204
                 << fion[is][i+2]
205
                 << " </force>\n  </atom>" << endl;
206

207 208 209
            i += 3;
          }
        }
210
        cout << "</atomset>" << endl;
211
      }
Francois Gygi committed
212
#if 1
Francois Gygi committed
213 214
      if ( s_.constraints.size() > 0 )
      {
Francois Gygi committed
215
        s_.constraints.compute_forces(mdionic_stepper->r0(), fion);
Francois Gygi committed
216 217
        if ( onpe0 )
        {
Francois Gygi committed
218
          s_.constraints.list_constraints(cout);
Francois Gygi committed
219 220 221
        }
      }
#endif
222
    }
223

Francois Gygi committed
224
    if ( onpe0 )
225 226
    {
      cout << "  <ekin_e> " << ekin_e << " </ekin_e>\n";
227

228 229 230 231 232 233 234 235 236
      if ( compute_forces )
      {
        cout << "  <ekin_ion> " << ekin_ion << " </ekin_ion>\n";
        cout << "  <temp_ion> " << mdionic_stepper->temp() << " </temp_ion>\n";
        cout << "  <eta_ion> " << mdionic_stepper->eta() << " </eta_ion>\n";
      }
      cout << "  <econst> " << energy+ekin_ion+ekin_e << " </econst>\n";
      cout << "  <ekin_ec> " << energy+ekin_ion+2*ekin_e << " </ekin_ec>\n";
    }
237

238 239 240 241
    if ( compute_stress )
    {
      compute_sigma();
      print_stress();
242

243 244 245
      if ( cell_dyn != "LOCKED" )
      {
        cell_stepper->compute_new_cell(sigma);
246

247 248
        // Update cell
        cell_stepper->update_cell();
249

250 251 252 253
        ef_.cell_moved();
        ef_.atoms_moved();
      }
    }
254

255 256 257 258
    if ( compute_forces )
    {
      ef_.atoms_moved();
    }
259

Francois Gygi committed
260
    tmap["charge"].start();
Francois Gygi committed
261
    cd_.update_density();
Francois Gygi committed
262
    tmap["charge"].stop();
Francois Gygi committed
263
    ef_.update_vhxc();
264
    energy =
265
      ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
266 267

    if ( s_.ctxt_.mype() == 0 )
268
      cout << "</iteration>" << endl;
269 270 271 272 273 274 275 276 277 278

    // print iteration time
    tm_iter.stop();
    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);
    if ( s_.ctxt_.myproc()==0 )
    {
279 280 281 282
      cout << "  <timing name=\"iteration\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
283
    }
Francois Gygi committed
284 285
    if ( compute_forces )
      s_.constraints.update_constraints(dt);
286
  } // iter
287

288 289
  // dwf and fion now contain the forces on wavefunctions and ions at the
  // endpoint
290

291
  mdwf_stepper->compute_wfv(dwf); // replace wfm by wfv
292

293 294
  if ( compute_forces )
  {
295
    // Note: next line function call updates velocities in the AtomSet
Francois Gygi committed
296
    mdionic_stepper->compute_v(energy,fion);
297
  }
298

299
  if ( cell_stepper != 0 ) delete cell_stepper;
300
}