CPSampleStepper.C 7.92 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
      string s = "name=\"" + (*i).first + "\"";
      cout << "<timing " << left << setw(22) << s
74 75
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
76
           << 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

138 139
  mdwf_stepper->compute_wfm(dwf);

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

147
    mdwf_stepper->update(dwf);
148

149
    ekin_e = mdwf_stepper->ekin();
150

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

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

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

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

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

209
    }
210

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

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

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

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

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

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

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

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

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

    // 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 )
    {
270 271
      string s = "name=\"iteration\"";
      cout << "<timing " << left << setw(22) << s
272 273
           << " min=\"" << setprecision(3) << tmin << "\""
           << " max=\"" << setprecision(3) << tmax << "\"/>"
274
           << 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
}