CPSampleStepper.C 8.31 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
3
// CPSampleStepper.C
4 5
//
////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
6
// $Id: CPSampleStepper.C,v 1.11 2005-09-16 23:08:11 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;

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
21 22
CPSampleStepper::CPSampleStepper(Sample& s) : 
  SampleStepper(s), cd_(s.wf), ef_(s,cd_), dwf(s.wf), wfv(s.wfv)
23 24 25 26 27 28 29 30 31
{
  mdwf_stepper = new MDWavefunctionStepper(s,tmap);
  assert(mdwf_stepper!=0);
  mdionic_stepper = 0;
  if ( s.ctrl.atoms_dyn != "LOCKED" )
    mdionic_stepper = new MDIonicStepper(s);
}

////////////////////////////////////////////////////////////////////////////////
32 33 34 35
CPSampleStepper::~CPSampleStepper(void)
{
  delete mdwf_stepper;
  if ( mdionic_stepper != 0 ) delete mdionic_stepper;
Francois Gygi committed
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
  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 )
    {
      cout << "<!-- timing "
           << setw(15) << (*i).first
           << " : " << setprecision(3) << setw(9) << tmin
           << " "   << setprecision(3) << setw(9) << tmax << " -->" << endl;
    }
  }
51 52 53 54
}

////////////////////////////////////////////////////////////////////////////////
void CPSampleStepper::step(int niter)
55
{
Francois Gygi committed
56
  const bool onpe0 = s_.ctxt_.onpe0();
57 58 59 60 61 62 63 64 65 66 67 68 69
  // 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;
  }
  
70 71 72 73 74 75 76 77 78
  AtomSet& atoms = s_.atoms;
  Wavefunction& wf = s_.wf;
  
  const double dt = s_.ctrl.dt;
  double ekin_ion=0.0,ekin_e, temp_ion, eta;
  
  const string wf_dyn = s_.ctrl.wf_dyn;
  assert(wf_dyn=="MD");
  const string atoms_dyn = s_.ctrl.atoms_dyn;
79
  const string cell_dyn = s_.ctrl.cell_dyn;
80 81 82
  
  const bool compute_hpsi = true;
  const bool compute_forces = ( atoms_dyn != "LOCKED" );
83
  const bool compute_stress = ( s_.ctrl.stress == "ON" );
84
  const bool use_confinement = ( s_.ctrl.ecuts > 0.0 );
85

86 87 88 89 90 91 92 93 94 95
  CellStepper* cell_stepper = 0;
  if ( cell_dyn == "SD" )
    cell_stepper = new SDCellStepper(s_);
    
  if ( s_.wfv == 0 )
  {
    s_.wfv = new Wavefunction(wf);
    s_.wfv->clear();
  }
  
Francois Gygi committed
96 97 98
  if ( mdionic_stepper )
    mdionic_stepper->setup_constraints();
  
99 100
  Timer tm_iter;
  
Francois Gygi committed
101
  tmap["charge"].start();
Francois Gygi committed
102
  cd_.update_density();
Francois Gygi committed
103 104
  tmap["charge"].stop();
  
Francois Gygi committed
105
  ef_.update_vhxc();
106
  double energy =
107
    ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
108
 
109 110
  mdwf_stepper->compute_wfm(dwf);

111 112 113 114 115 116
  for ( int iter = 0; iter < niter; iter++ )
  {
    tm_iter.reset();
    tm_iter.start();
    if ( s_.ctxt_.mype() == 0 )
      cout << "  <iteration count=\"" << iter+1 << "\">\n";
Francois Gygi committed
117 118 119
      
    if ( mdionic_stepper )
      atoms.sync();
120
 
121
    mdwf_stepper->update(dwf);
122 123 124
      
    ekin_e = mdwf_stepper->ekin();
 
Francois Gygi committed
125
    if ( onpe0 )
126 127 128
    {
      cout.setf(ios::fixed,ios::floatfield);
      cout.setf(ios::right,ios::adjustfield);
129
      cout << "  <ekin>     " << setprecision(8)
130
           << setw(15) << ef_.ekin() << " </ekin>\n";
131
      if ( use_confinement )
132
      {
133 134
        cout << "  <econf>    " << setw(15) << ef_.econf()
             << " </econf>\n";
135
      }
136 137 138 139 140 141 142
      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"
143
           << flush;
144 145 146 147 148 149 150 151 152
      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;
      }
153 154 155 156
    }
 
    if ( compute_forces )
    {
157 158
      if ( iter > 0 ) 
      {
Francois Gygi committed
159
        mdionic_stepper->compute_v(energy,fion);
160
      }
Francois Gygi committed
161
      mdionic_stepper->compute_r(energy,fion);
162 163
      ekin_ion = mdionic_stepper->ekin();
      
Francois Gygi committed
164
      if ( onpe0 )
165 166 167 168 169 170 171 172 173 174 175
      {
        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> " 
176 177 178 179 180 181 182 183 184 185 186
                 << mdionic_stepper->r0(is,i) << " "
                 << mdionic_stepper->r0(is,i+1) << " " 
                 << mdionic_stepper->r0(is,i+2) << " </position>\n"
                 << "    <velocity> " 
                 << mdionic_stepper->v0(is,i) << " "
                 << mdionic_stepper->v0(is,i+1) << " " 
                 << mdionic_stepper->v0(is,i+2) << " </velocity>\n"
                 << "    <force> " 
                 << fion[is][i] << " "
                 << fion[is][i+1] << " " 
                 << fion[is][i+2]
187 188 189 190 191 192
                 << " </force>\n  </atom>" << endl;
 
            i += 3;
          }
        }
      }
Francois Gygi committed
193
#if 1
Francois Gygi committed
194 195
      if ( s_.constraints.size() > 0 )
      {
Francois Gygi committed
196
        s_.constraints.compute_forces(mdionic_stepper->r0(), fion);
Francois Gygi committed
197 198
        if ( onpe0 )
        {
Francois Gygi committed
199
          s_.constraints.list_constraints(cout);
Francois Gygi committed
200 201 202
        }
      }
#endif
203 204
    }
 
Francois Gygi committed
205
    if ( onpe0 )
206 207 208 209 210 211 212 213 214 215 216 217 218
    {
      cout << "  <ekin_e> " << ekin_e << " </ekin_e>\n";
 
      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";
    }
    
219 220
    if ( compute_stress )
    {
Francois Gygi committed
221
      if ( onpe0 )            
222 223 224 225 226
      {                                  
        cout << "<unit_cell>" << endl;   
        cout << s_.wf.cell();            
        cout << "</unit_cell>" << endl;  
      }                                  
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
      compute_sigma();
      print_stress();
      
      if ( cell_dyn != "LOCKED" )
      {
        cell_stepper->compute_new_cell(sigma);
      
        // Update cell
        cell_stepper->update_cell();
      
        ef_.cell_moved();
        ef_.atoms_moved();
      }
    }
    
    if ( compute_forces )
    {
      ef_.atoms_moved();
    }
    
Francois Gygi committed
247
    tmap["charge"].start();
Francois Gygi committed
248
    cd_.update_density();
Francois Gygi committed
249
    tmap["charge"].stop();
Francois Gygi committed
250
    ef_.update_vhxc();
251
    energy =
252
      ef_.energy(compute_hpsi,dwf,compute_forces,fion,compute_stress,sigma_eks);
253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270

    if ( s_.ctxt_.mype() == 0 )
      cout << "  </iteration>" << endl;

    // 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 )
    {
      cout << "  <!-- timing "
           << setw(15) << "iteration"
           << " : " << setprecision(3) << setw(9) << tmin
           << " "   << setprecision(3) << setw(9) << tmax << " -->" << endl;
    }
Francois Gygi committed
271 272
    if ( compute_forces )
      s_.constraints.update_constraints(dt);
273 274 275 276 277
  } // iter
 
  // dwf and fion now contain the forces on wavefunctions and ions at the
  // endpoint
 
278
  mdwf_stepper->compute_wfv(dwf); // replace wfm by wfv
279 280 281
  
  if ( compute_forces )
  {
282
    // Note: next line function call updates velocities in the AtomSet
Francois Gygi committed
283
    mdionic_stepper->compute_v(energy,fion);
284
  }
285 286
  
  if ( cell_stepper != 0 ) delete cell_stepper;
287
}