CPSampleStepper.C 7.95 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
// CPSampleStepper.C
16 17 18 19 20 21 22
//
////////////////////////////////////////////////////////////////////////////////

#include "CPSampleStepper.h"
#include "SlaterDet.h"
#include "MDWavefunctionStepper.h"
#include "MDIonicStepper.h"
23
#include "SDCellStepper.h"
24
#include "CGCellStepper.h"
25
#include "Basis.h"
Francois Gygi committed
26
#include "Species.h"
27 28 29 30 31 32

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

////////////////////////////////////////////////////////////////////////////////
33
CPSampleStepper::CPSampleStepper(Sample& s) :
Francois Gygi committed
34
  SampleStepper(s), cd_(s.wf), ef_(s,cd_), dwf(s.wf), wfv(s.wfv)
35
{
36 37
  const double emass = s.ctrl.emass;
  const double dt = s.ctrl.dt;
38 39 40 41 42 43 44 45
  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;
  }
46 47 48 49 50
  if ( s.wfv == 0 )
  {
    s.wfv = new Wavefunction(s.wf);
    s.wfv->clear();
  }
51
  mdwf_stepper = new MDWavefunctionStepper(s.wf,s.wfv,dt,dt2bye,tmap);
52 53 54 55 56 57 58
  assert(mdwf_stepper!=0);
  mdionic_stepper = 0;
  if ( s.ctrl.atoms_dyn != "LOCKED" )
    mdionic_stepper = new MDIonicStepper(s);
}

////////////////////////////////////////////////////////////////////////////////
59 60 61 62
CPSampleStepper::~CPSampleStepper(void)
{
  delete mdwf_stepper;
  if ( mdionic_stepper != 0 ) delete mdionic_stepper;
Francois Gygi committed
63 64 65 66 67 68 69 70 71
  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 )
    {
72 73 74 75 76
      cout << "<timing name=\""
           << setw(15) << (*i).first << "\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
Francois Gygi committed
77 78
    }
  }
79 80 81 82
}

////////////////////////////////////////////////////////////////////////////////
void CPSampleStepper::step(int niter)
83
{
Francois Gygi committed
84
  const bool onpe0 = s_.ctxt_.onpe0();
85 86 87 88 89 90 91 92 93 94 95 96
  // 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;
  }
97

98 99
  AtomSet& atoms = s_.atoms;
  Wavefunction& wf = s_.wf;
100

101
  const double dt = s_.ctrl.dt;
102
  double ekin_ion=0.0,ekin_e;
103

104 105 106
  const string wf_dyn = s_.ctrl.wf_dyn;
  assert(wf_dyn=="MD");
  const string atoms_dyn = s_.ctrl.atoms_dyn;
107
  const string cell_dyn = s_.ctrl.cell_dyn;
108

109 110
  const bool compute_hpsi = true;
  const bool compute_forces = ( atoms_dyn != "LOCKED" );
111
  const bool compute_stress = ( s_.ctrl.stress == "ON" );
112

113 114 115
  CellStepper* cell_stepper = 0;
  if ( cell_dyn == "SD" )
    cell_stepper = new SDCellStepper(s_);
116 117
  else if ( cell_dyn == "CG" )
    cell_stepper = new CGCellStepper(s_);
118

119 120 121 122 123
  if ( s_.wfv == 0 )
  {
    s_.wfv = new Wavefunction(wf);
    s_.wfv->clear();
  }
124

Francois Gygi committed
125 126
  if ( mdionic_stepper )
    mdionic_stepper->setup_constraints();
127

128
  Timer tm_iter;
129

Francois Gygi committed
130
  tmap["charge"].start();
Francois Gygi committed
131
  cd_.update_density();
Francois Gygi committed
132
  tmap["charge"].stop();
133

134
  ef_.update_vhxc(compute_stress);
135
  double energy =
136
    ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
137
  double enthalpy = ef_.enthalpy();
138

139 140
  mdwf_stepper->compute_wfm(dwf);

141 142 143 144 145
  for ( int iter = 0; iter < niter; iter++ )
  {
    tm_iter.reset();
    tm_iter.start();
    if ( s_.ctxt_.mype() == 0 )
146
      cout << "<iteration count=\"" << iter+1 << "\">\n";
147

148
    mdwf_stepper->update(dwf);
149

150
    ekin_e = mdwf_stepper->ekin();
151

Francois Gygi committed
152
    if ( onpe0 )
153
    {
154 155 156
      cout << ef_;
      if ( ef_.el_enth() )
        cout << *ef_.el_enth();
157
    }
158

159 160
    if ( compute_forces )
    {
161
      if ( iter > 0 )
162
      {
Francois Gygi committed
163
        mdionic_stepper->compute_v(energy,fion);
164
      }
Francois Gygi committed
165
      mdionic_stepper->compute_r(energy,fion);
166
      ekin_ion = mdionic_stepper->ekin();
167

Francois Gygi committed
168
      if ( onpe0 )
169
      {
170
        cout << "<atomset>" << endl;
Francois Gygi committed
171
        cout << atoms.cell();
172 173 174 175 176 177 178 179 180
        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"
181
                 << "    <position> "
182
                 << mdionic_stepper->r0(is,i) << " "
183
                 << mdionic_stepper->r0(is,i+1) << " "
184
                 << mdionic_stepper->r0(is,i+2) << " </position>\n"
185
                 << "    <velocity> "
186
                 << mdionic_stepper->v0(is,i) << " "
187
                 << mdionic_stepper->v0(is,i+1) << " "
188
                 << mdionic_stepper->v0(is,i+2) << " </velocity>\n"
189
                 << "    <force> "
190
                 << fion[is][i] << " "
191
                 << fion[is][i+1] << " "
192
                 << fion[is][i+2]
193
                 << " </force>\n  </atom>" << endl;
194

195 196 197
            i += 3;
          }
        }
198
        cout << "</atomset>" << endl;
199
      }
200

Francois Gygi committed
201 202
      if ( s_.constraints.size() > 0 )
      {
Francois Gygi committed
203
        s_.constraints.compute_forces(mdionic_stepper->r0(), fion);
Francois Gygi committed
204 205
        if ( onpe0 )
        {
Francois Gygi committed
206
          s_.constraints.list_constraints(cout);
Francois Gygi committed
207 208
        }
      }
209

210
    }
211

Francois Gygi committed
212
    if ( onpe0 )
213 214
    {
      cout << "  <ekin_e> " << ekin_e << " </ekin_e>\n";
215

216 217 218 219 220 221
      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";
      }
222 223 224
      double econst = energy + ekin_ion + ekin_e;
      if ( mdionic_stepper )
        econst += mdionic_stepper->ekin_stepper();
225 226
      cout << "  <econst> " << econst << " </econst>\n";
      cout << "  <ekin_ec> " << econst + ekin_e << " </ekin_ec>\n";
227
    }
228

229 230 231 232
    if ( compute_stress )
    {
      compute_sigma();
      print_stress();
233

234 235
      if ( cell_dyn != "LOCKED" )
      {
236
        cell_stepper->compute_new_cell(energy,sigma,fion);
237

238 239
        // Update cell
        cell_stepper->update_cell();
240

241 242 243 244
        ef_.cell_moved();
        ef_.atoms_moved();
      }
    }
245

246 247 248 249
    if ( compute_forces )
    {
      ef_.atoms_moved();
    }
250

Francois Gygi committed
251
    tmap["charge"].start();
Francois Gygi committed
252
    cd_.update_density();
Francois Gygi committed
253
    tmap["charge"].stop();
254
    ef_.update_vhxc(compute_stress);
255
    energy =
256
      ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
257
    enthalpy = ef_.enthalpy();
258 259

    if ( s_.ctxt_.mype() == 0 )
260
      cout << "</iteration>" << endl;
261 262 263 264 265 266 267 268 269 270

    // 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 )
    {
271 272 273 274
      cout << "  <timing name=\"iteration\""
           << " min=\"" << setprecision(3) << setw(9) << tmin << "\""
           << " max=\"" << setprecision(3) << setw(9) << tmax << "\"/>"
           << endl;
275
    }
Francois Gygi committed
276 277
    if ( compute_forces )
      s_.constraints.update_constraints(dt);
278
  } // iter
279

280 281
  // dwf and fion now contain the forces on wavefunctions and ions at the
  // endpoint
282

283
  mdwf_stepper->compute_wfv(dwf); // replace wfm by wfv
284

285 286
  if ( compute_forces )
  {
287
    // Note: next line function call updates velocities in the AtomSet
Francois Gygi committed
288
    mdionic_stepper->compute_v(energy,fion);
289
  }
290

291
  if ( cell_stepper != 0 ) delete cell_stepper;
292
}