CGCellStepper.C 5.52 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2011 The Regents of the University of California
//
// This file is part of Qbox
//
// Qbox is distributed under the terms of the GNU General Public License
// as published by the Free Software Foundation, either version 2 of
// 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// CGCellStepper.C
//
////////////////////////////////////////////////////////////////////////////////

#include "CGCellStepper.h"
20
#include "CGOptimizer.h"
21 22 23
using namespace std;

////////////////////////////////////////////////////////////////////////////////
Francois Gygi committed
24
CGCellStepper::CGCellStepper(Sample& s) : CellStepper(s),
25
  cgopt_(CGOptimizer(3*s.atoms.size()+9)), cell0(s_.atoms.cell())
26 27
{
  nat_ = atoms_.size();
28 29 30
  cgopt_.set_alpha_start(0.05);
  cgopt_.set_alpha_max(1.0);
  cgopt_.set_beta_max(10.0);
31
#ifdef DEBUG
32 33
  if ( s.ctxt_.onpe0() )
    cgopt_.set_debug_print();
34
#endif
35

36 37 38
  rp_.resize(s.atoms.nsp());
  for ( int is = 0; is < rp_.size(); is++ )
    rp_[is].resize(3*atoms_.na(is));
39

40 41 42 43
  // store full strain tensor u for consistency of dot products
  u_.resize(9);
  up_.resize(9);
  for ( int i = 0; i < u_.size(); i++ )
44
  {
45 46
    u_[i] = 0.0;
    up_[i] = 0.0;
47 48 49 50
  }
}

////////////////////////////////////////////////////////////////////////////////
51
void CGCellStepper::compute_new_cell(double e, const valarray<double>& sigma,
52 53 54 55
  const std::vector<std::vector< double> >& fion)
{
  // compute new cell and ionic positions using the stress tensor sigma
  // and the forces on ions fion
56
  const UnitCell cell = s_.wf.cell();
57 58

  // total number of dofs: 3* natoms + cell parameters
59
  valarray<double> x(3*nat_+9), xp(3*nat_+9), g(3*nat_+9);
60

61
  // copy current positions into x
62
  vector<vector<double> > r0, gvec;
63 64 65
  atoms_.get_positions(r0);
  double tmp3[3];
  // convert position from r0 to tau (relative) coordinates: tau = A^-1 R
66
  for ( int is = 0, i = 0; is < r0.size(); is++ )
67 68 69 70
  {
    for ( int ia = 0; ia < r0[is].size()/3; ia++ )
    {
      cell.vecmult3x3(cell.amat_inv(),&r0[is][3*ia],tmp3);
71 72 73
      x[i++]=tmp3[0];
      x[i++]=tmp3[1];
      x[i++]=tmp3[2];
74 75 76
    }
  }

77 78 79 80
  // copy current strain tensor into x
  for ( int i = 0; i < 9; i++ )
    x[3*nat_+i] = u_[i];

81
  // convert forces on positions to forces on tau coordinates, store -f in gvec
Francois Gygi committed
82
  // f = A^-1 * fion
83 84
  gvec.resize(r0.size());
  for ( int is = 0; is < r0.size(); is++ )
85
  {
86 87
    gvec[is].resize(r0[is].size());
    for ( int ia = 0; ia < r0[is].size()/3; ia++ )
88
    {
Francois Gygi committed
89
      cell.vecmult3x3(cell.amat_inv(),&fion[is][3*ia],tmp3);
90 91 92
      gvec[is][3*ia+0]=-tmp3[0];
      gvec[is][3*ia+1]=-tmp3[1];
      gvec[is][3*ia+2]=-tmp3[2];
93 94 95
    }
  }

96 97 98 99 100 101 102
  // project the gradient gvec in a direction compatible with constraints
  s_.constraints.enforce_v(r0,gvec);

  for ( int j = 0, is = 0; is < r0.size(); is++ )
    for ( int i = 0; i < r0[is].size(); i++ )
      g[j++] = gvec[is][i];

103 104 105 106 107
  // compute descent direction in strain space
  // gradient g = - sigma * volume
  g[3*nat_+0] = -sigma[0] * cell.volume();
  g[3*nat_+1] = -sigma[3] * cell.volume();
  g[3*nat_+2] = -sigma[5] * cell.volume();
108

109 110 111
  g[3*nat_+3] = -sigma[3] * cell.volume();
  g[3*nat_+4] = -sigma[1] * cell.volume();
  g[3*nat_+5] = -sigma[4] * cell.volume();
112

113 114 115
  g[3*nat_+6] = -sigma[5] * cell.volume();
  g[3*nat_+7] = -sigma[4] * cell.volume();
  g[3*nat_+8] = -sigma[2] * cell.volume();
116

117
  // the vector g now contains the gradient of the energy in tau+strain space
118

119 120
  // enforce constraints on the unit cell
  // Next line only affects the cell components of g
121
  enforce_constraints(&g[3*nat_]);
Francois Gygi committed
122

123 124 125
  // the vector g now contains a gradient of the energy in tau+strain space
  // compatible with atoms and cell constraints

126
  // CG algorithm
127

128
  cgopt_.compute_xp(x,e,g,xp);
Francois Gygi committed
129 130 131

  if ( s_.ctxt_.onpe0() )
  {
132
    cout << "  CGCellStepper: alpha = " << cgopt_.alpha() << endl;
Francois Gygi committed
133
  }
134

Francois Gygi committed
135
  for ( int i = 0; i < 9; i++ )
Francois Gygi committed
136
    up_[i] = xp[3*nat_+i];
137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158

  // enforce cell_lock constraints
  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);
Francois Gygi committed
159

Francois Gygi committed
160
  // compute new atomic positions rp_[is][3*ia+j] from xp and cellp
161 162 163 164 165 166 167
  for ( int is = 0, i = 0; is < fion.size(); is++ )
  {
    for ( int ia = 0; ia < fion[is].size()/3; ia++ )
    {
      tmp3[0] = xp[i+0];
      tmp3[1] = xp[i+1];
      tmp3[2] = xp[i+2];
168
      cellp.vecmult3x3(cellp.amat(),tmp3,&rp_[is][3*ia]);
169 170 171 172 173
      i+=3;
    }
  }

  // enforce constraints on atomic positions
174
  s_.constraints.enforce_r(r0,rp_);
175 176 177 178 179
}

////////////////////////////////////////////////////////////////////////////////
void CGCellStepper::update_cell(void)
{
180
  s_.atoms.sync_positions(rp_);
181
  s_.atoms.set_positions(rp_);
182
  s_.atoms.sync_cell(cellp);
183
  s_.atoms.set_cell(cellp);
184
  u_ = up_;
185 186 187 188 189 190

  // 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());
}