SDCellStepper.C 3.87 KB
Newer Older
1 2
////////////////////////////////////////////////////////////////////////////////
//
3
// Copyright (c) 2011 The Regents of the University of California
Francois Gygi committed
4 5 6
//
// 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
// SDCellStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "SDCellStepper.h"
20
using namespace std;
21 22

////////////////////////////////////////////////////////////////////////////////
23
SDCellStepper::SDCellStepper(Sample& s) : CellStepper(s)
24
{
25 26 27 28 29 30 31 32 33
  // strain tensors u, up
  u_.resize(9);
  up_.resize(9);
  for ( int i = 0; i < u_.size(); i++ )
  {
    u_[i] = 0.0;
    up_[i] = 0.0;
  }
}
34

35 36 37 38 39 40
////////////////////////////////////////////////////////////////////////////////
void SDCellStepper::compute_new_cell(double e0, const valarray<double>& sigma,
  const std::vector<std::vector< double> >& fion)
{
  // compute new cell and ionic positions using the stress tensor sigma
  const UnitCell cell0 = s_.wf.cell();
41
  const double cell_mass = s_.ctrl.cell_mass;
42

43 44
  if ( cell_mass <= 0.0 )
  {
45
    cellp = cell0;
46 47
    if ( s_.ctxt_.onpe0() )
    {
48 49
      cout << " SDCellStepper::compute_new_cell: cell mass is zero\n"
           << "     cannot update cell" << endl;
50 51 52
      return;
    }
  }
53

54
  assert(sigma.size()==6);
55 56
  const double dt = s_.ctrl.dt;
  const double dt2bym = dt*dt/cell_mass;
57

58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
  double g[9];
  // compute descent direction in strain space
  const double omega = cell0.volume();
  // gradient g = - sigma * volume
  g[0] = -sigma[0] * omega;
  g[1] = -sigma[3] * omega;
  g[2] = -sigma[5] * omega;

  g[3] = -sigma[3] * omega;
  g[4] = -sigma[1] * omega;
  g[5] = -sigma[4] * omega;

  g[6] = -sigma[5] * omega;
  g[7] = -sigma[4] * omega;
  g[8] = -sigma[2] * omega;

  // the vector g now contains the gradient of the energy in strain space

  // SD algorithm
  for ( int i = 0; i < 9; i++ )
    up_[i] = - dt2bym * g[i];
Francois Gygi committed
79

80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
  // check for cell_lock constraints and modify up_ if needed
  enforce_constraints(&up_[0]);

  // compute cellp = ( I + Up ) * A0
  // symmetric matrix I+U stored in iumat: xx, yy, zz, xy, yz, xz
  double iupmat[6];
  iupmat[0] = 1.0 + up_[0];
  iupmat[1] = 1.0 + up_[4];
  iupmat[2] = 1.0 + up_[8];
  iupmat[3] = 0.5 * ( up_[1] + up_[3] );
  iupmat[4] = 0.5 * ( up_[5] + up_[7] );
  iupmat[5] = 0.5 * ( up_[2] + up_[6] );

  const double *a0mat = cell0.amat();
  double apmat[9];
  cell0.smatmult3x3(&iupmat[0],a0mat,apmat);

  D3vector a0p(apmat[0],apmat[1],apmat[2]);
  D3vector a1p(apmat[3],apmat[4],apmat[5]);
  D3vector a2p(apmat[6],apmat[7],apmat[8]);
  cellp.set(a0p,a1p,a2p);
101 102 103 104 105 106
}

////////////////////////////////////////////////////////////////////////////////
void SDCellStepper::update_cell(void)
{
  const UnitCell& cell = s_.wf.cell();
107

108
  // rescale atomic positions in AtomSet
109

110 111 112
  // r_new = A_new A_old^-1 r_old
  vector<vector<double> > r;
  s_.atoms.get_positions(r);
113

114 115 116 117 118 119 120 121 122 123 124 125 126
  double tau[3];
  for ( int is = 0; is < r.size(); is++ )
  {
    // transform r to tau: multiply by A^-1
    const int nais = r[is].size()/3;
    for ( int ia = 0; ia < nais; ia++ )
    {
      // multiply r[is][ia] by A_old^-1, result in tau
      cell.vecmult3x3(cell.amat_inv(),&r[is][3*ia],&tau[0]);
      // multiply tau by A_new, result in r[is][3*ia]
      cellp.vecmult3x3(cellp.amat(),&tau[0],&r[is][3*ia]);
    }
  }
127
  s_.atoms.sync_positions(r);
128
  s_.atoms.set_positions(r);
129
  s_.atoms.sync_cell(cellp);
Francois Gygi committed
130
  s_.atoms.set_cell(cellp);
131

132 133 134 135 136
  // resize wavefunction and basis sets
  s_.wf.resize(cellp,s_.wf.refcell(),s_.wf.ecut());
  if ( s_.wfv != 0 )
    s_.wfv->resize(cellp,s_.wf.refcell(),s_.wf.ecut());
}