SDCellStepper.C 4.94 KB
Newer Older
1 2 3 4 5
////////////////////////////////////////////////////////////////////////////////
//
// SDCellStepper.C
//
////////////////////////////////////////////////////////////////////////////////
6
// $Id: SDCellStepper.C,v 1.3 2004-05-04 21:24:11 fgygi Exp $
7 8 9 10 11 12

#include "SDCellStepper.h"

////////////////////////////////////////////////////////////////////////////////
void SDCellStepper::compute_new_cell(const valarray<double>& sigma)
{
13
  //cout << " SDCellStepper::compute_new_cell" << endl;
14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
  // multiply stress by A^T to get dE/da_ij
  valarray<double> deda(9);
  
  const UnitCell& cell = s_.wf.cell();
  const double cell_mass = s_.ctrl.cell_mass;
  
  if ( cell_mass <= 0.0 )
  {
    if ( s_.ctxt_.onpe0() )
    {
      cout << "<!-- SDCellStepper::compute_new_cell: cell mass is zero\n"
           << "     cannot update cell -->" << endl;
      return;
    }
  }
  
  // deda = - omega * sigma * A^-T
  cell.compute_deda(sigma,deda);
  
33
  //cout << " SDCellStepper: cell derivatives before constraints" << endl;
34 35
  //for ( int i = 0; i < 9; i++ )
  //  cout << " deda[" << i << "] = " << deda[i] << endl;
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
  
  string cell_lock = s_.ctrl.cell_lock;
  if ( cell_lock != "OFF" )
  {
    // constraints on the cell derivatives
    if ( cell_lock.find("A") != string::npos )
    {
      // vector A is locked
      deda[0] = deda[1] = deda[2] = 0.0;
    }
    if ( cell_lock.find("B") != string::npos )
    {
      // vector B is locked
      deda[3] = deda[4] = deda[5] = 0.0;
    }
    if ( cell_lock.find("C") != string::npos )
    {
      // vector C is locked
      deda[6] = deda[7] = deda[8] = 0.0;
    }
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
    
    // Check is cell shape should be preserved (if "S" is present in cell_lock)
    // The only changes allowed are renormalizations of a,b,c
    if ( cell_lock.find("S") != string::npos )
    {
      // projection of d in the direction of e
      D3vector d,e;
      
      d = D3vector(deda[0],deda[1],deda[2]);
      e = cell.a(0) / length(cell.a(0));
      d = (d * e) * e;
      deda[0] = d.x; deda[1] = d.y; deda[2] = d.z;
      
      d = D3vector(deda[3],deda[4],deda[5]);
      e = cell.a(1) / length(cell.a(1));
      d = (d * e) * e;
      deda[3] = d.x; deda[4] = d.y; deda[5] = d.z;
      
      d = D3vector(deda[6],deda[7],deda[8]);
      e = cell.a(2) / length(cell.a(2));
      d = (d * e) * e;
      deda[6] = d.x; deda[7] = d.y; deda[8] = d.z;
    }
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
    
    if ( cell_lock == "R" )
    {
      // preserve aspect ratio
      // deda must be proportional to A
      // All vectors are rescaled by the same constant
      // rescale cell by coefficient alpha, i.e. 
      // deda = alpha * A
      // where alpha * A is the projection of dE/dA in the direction of A
      // alpha = tr (A^T * deda) / || A ||^2  (matrix scalar product)
      const double *a = cell.amat();
      const double num = a[0]*deda[0] + a[1]*deda[1] + a[2]*deda[2] +
                         a[3]*deda[3] + a[4]*deda[4] + a[5]*deda[5] +
                         a[6]*deda[6] + a[7]*deda[7] + a[8]*deda[8];
      const double denom = a[0]*a[0] + a[1]*a[1] + a[2]*a[2] +
                           a[3]*a[3] + a[4]*a[4] + a[5]*a[5] +
                           a[6]*a[6] + a[7]*a[7] + a[8]*a[8];
      const double alpha = num / denom;
      deda = valarray<double>(a,9);
      deda *= alpha;
    }
    //cout << " SDCellStepper: cell derivatives after constraints" << endl;
101 102
    //for ( int i = 0; i < 9; i++ )
    //  cout << " deda[" << i << "] = " << deda[i] << endl;
103 104 105 106 107 108 109 110 111 112 113 114
  }
  
  const double dt = s_.ctrl.dt;
  const double dt2bym = dt*dt/cell_mass;
  
  // cellp = cell - deda * dt^2 / cell_mass
  D3vector a0p = cell.a(0) - dt2bym * D3vector(deda[0],deda[1],deda[2]);
  D3vector a1p = cell.a(1) - dt2bym * D3vector(deda[3],deda[4],deda[5]);
  D3vector a2p = cell.a(2) - dt2bym * D3vector(deda[6],deda[7],deda[8]);
  
  cellp = UnitCell(a0p,a1p,a2p);
  
115 116
  //cout << " SDCellStepper::compute_new_cell: cellp: " << endl;
  //cout << cellp;
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
  
}

////////////////////////////////////////////////////////////////////////////////
void SDCellStepper::update_cell(void)
{
  const UnitCell& cell = s_.wf.cell();
  
  // rescale atomic positions in AtomSet
  
  // r_new = A_new A_old^-1 r_old
  vector<vector<double> > r;
  s_.atoms.get_positions(r);
  const double* const ainv = cell.amat_inv();
  const double* const ap = cellp.amat();
  
  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]);
    }
  }
  s_.atoms.set_positions(r);
  
  // resize wavefunction and basis sets
  
150
  //cout << " SDCellStepper::update_cell" << endl;
151 152 153 154 155 156 157
  s_.wf.resize(cellp,s_.wf.refcell(),s_.wf.ecut());
  if ( s_.wfv != 0 )
  {
    s_.wfv->resize(cellp,s_.wf.refcell(),s_.wf.ecut());
    // s_.wfv->clear();
  }
}