StrainCmd.h 4.97 KB
Newer Older
Francois Gygi committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
////////////////////////////////////////////////////////////////////////////////
//
// Copyright (c) 2008 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/>.
//
////////////////////////////////////////////////////////////////////////////////
//
// StrainCmd.h:
//
////////////////////////////////////////////////////////////////////////////////

#ifndef STRAINCMD_H
#define STRAINCMD_H

#include <string>
#include <cstdlib>
#include <iostream>
using namespace std;

#include "UserInterface.h"
#include "Sample.h"

class StrainCmd : public Cmd
{
  public:

  Sample *s;

  StrainCmd(Sample *sample) : s(sample) {};

38 39
  const char *name(void) const { return "strain"; }
  const char *help_msg(void) const
Francois Gygi committed
40 41 42
  {
    return
    "\n strain\n\n"
Francois Gygi committed
43
    " syntax: strain [-atomsonly] [-inverse] uxx uyy uzz uxy uyz uxz \n\n"
Francois Gygi committed
44 45 46 47 48 49 50 51 52 53
    "   The strain command deforms the unit cell according to\n"
    "   the given strain tensor u. Atomic positions are modified\n"
    "   correspondingly. If -atomsonly is specified, only atomic\n"
    "   positions are modified and the cell is unchanged. If the\n"
    "   -inverse flag is used, the inverse transformation is applied.\n\n";
  }

  int action(int argc, char **argv)
  {
    const string usage =
Francois Gygi committed
54
    " use: strain [-atomsonly] [-inverse] uxx uyy uzz uxy uyz uxz";
Francois Gygi committed
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 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 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189
    // strain must have 7, 8 or 9 arguments including the command name
    if ( argc < 7 || argc > 9 )
    {
      if ( ui->onpe0() )
        cout << usage << endl;
      return 1;
    }

    bool atomsonly = false;
    bool inverse = false;
    int iu = 1;
    for ( int iarg = 0; iarg < 2; iarg++ )
    {
      string arg = argv[iarg];
      if ( arg == "-atomsonly" )
      {
        atomsonly = true;
        iu++;
      }
      if ( arg == "-inverse" )
      {
        inverse = true;
        iu++;
      }
    }
    if ( argc != iu+6 )
    {
      if ( ui->onpe0() )
        cout << usage << endl;
      return 1;
    }

    vector<double> u(6);
    for ( int i = 0; i < 6; i++ )
      u[i] = atof(argv[iu+i]);

#ifdef DEBUG
    if ( ui->onpe0() )
    {
      cout << " u: input: " << endl;
      cout << u[0] << " " << u[3] << " " << u[5] << endl;
      cout << u[3] << " " << u[1] << " " << u[4] << endl;
      cout << u[5] << " " << u[4] << " " << u[2] << endl;
    }
#endif

    if ( inverse )
    {
      // replace u by inv(I+u)
      const double a0 = u[0] + 1.0;
      const double a1 = u[1] + 1.0;
      const double a2 = u[2] + 1.0;
      const double a3 = u[3];
      const double a4 = u[4];
      const double a5 = u[5];

      const double det = a0*(a1*a2-a4*a4)-a3*(a3*a2-a5*a4)+a5*(a3*a4-a5*a1);
      if ( fabs(det) < 1.e-8 )
      {
        if ( ui->onpe0() )
        {
          cout << " transformation is near-singular: det(I+u) < 1.e-8" << endl;
        }
        return 1;
      }
      const double detinv = 1.0 / det;
      // replace u with (a^-1 - I)
      u[0] = detinv * ( a1*a2-a4*a4 ) - 1.0;
      u[1] = detinv * ( a0*a2-a5*a5 ) - 1.0;
      u[2] = detinv * ( a0*a1-a3*a3 ) - 1.0;
      u[3] = detinv * ( a4*a5-a3*a2 );
      u[4] = detinv * ( a5*a3-a0*a4 );
      u[5] = detinv * ( a3*a4-a5*a1 );

    }
#ifdef DEBUG
    if ( ui->onpe0() )
    {
      cout << " u: used: " << endl;
      cout << u[0] << " " << u[3] << " " << u[5] << endl;
      cout << u[3] << " " << u[1] << " " << u[4] << endl;
      cout << u[5] << " " << u[4] << " " << u[2] << endl;
    }
#endif

    UnitCell cell = s->atoms.cell();
    if ( !atomsonly )
    {
      // compute cell deformation: damat = u * amat
      vector<double> damat(9);

      cell.smatmult3x3(&u[0], cell.amat(), &damat[0]);

      D3vector a0new(cell.a(0)+D3vector(damat[0],damat[1],damat[2]));
      D3vector a1new(cell.a(1)+D3vector(damat[3],damat[4],damat[5]));
      D3vector a2new(cell.a(2)+D3vector(damat[6],damat[7],damat[8]));
      cell.set(a0new,a1new,a2new);

      s->atoms.set_cell(a0new,a1new,a2new);

      s->wf.resize(cell,s->wf.refcell(),s->wf.ecut());
      if ( s->wfv != 0 )
      {
        s->wfv->resize(cell,s->wf.refcell(),s->wf.ecut());
        s->wfv->clear();
      }

      if ( ui->onpe0() )
      {
        cout << s->atoms.cell();
      }
    }

    // atomic displacements
    vector<vector<double> > tau;
    s->atoms.get_positions(tau);

    vector<vector<double> > dtau;
    dtau.resize(tau.size());
    for ( int is = 0; is < dtau.size(); is++ )
      dtau[is].resize(tau[is].size());

    for ( int is = 0; is < tau.size(); is++ )
      for ( int ia = 0; ia < s->atoms.na(is); ia++ )
        cell.vecsmult3x3(&u[0], &tau[is][3*ia], &dtau[is][3*ia]);
    for ( int is = 0; is < tau.size(); is++ )
      for ( int i = 0; i < tau[is].size(); i++ )
        tau[is][i] += dtau[is][i];

    s->atoms.set_positions(tau);

    return 0;
  }
};
#endif