rel1_45_0

git-svn-id: http://qboxcode.org/svn/qb/trunk@659 cba15fb0-1239-40c8-b417-11db7ca47a34

src/Makefile

@@ -11,7 +11,7 @@
 # or <http://www.gnu.org/licenses/>.
 #
 #-------------------------------------------------------------------------------
-# $Id: Makefile,v 1.57 2008-09-08 15:56:18 fgygi Exp $
+# $Id: Makefile,v 1.58 2008-09-15 15:00:31 fgygi Exp $
 #------------------------------------------------------------------------------
 #
 include $(TARGET).mk
@@ -243,7 +243,8 @@ BOSampleStepper.o: PSDWavefunctionStepper.h PSDAWavefunctionStepper.h
 BOSampleStepper.o: SDIonicStepper.h IonicStepper.h Species.h
 BOSampleStepper.o: SDAIonicStepper.h LineMinimizer.h CGIonicStepper.h
 BOSampleStepper.o: MDIonicStepper.h SDCellStepper.h CellStepper.h
-BOSampleStepper.o: Preconditioner.h AndersonMixer.h
+BOSampleStepper.o: Preconditioner.h AndersonMixer.h MLWFTransform.h
+BOSampleStepper.o: BasisMapping.h
 BOSampleStepper.o: SampleStepper.h Timer.h EnergyFunctional.h ChargeDensity.h
 BOSampleStepper.o: Context.h StructureFactor.h Sample.h AtomSet.h Atom.h
 BOSampleStepper.o: D3vector.h UnitCell.h ConstraintSet.h Wavefunction.h
@@ -352,7 +353,7 @@ IonicStepper.o: Species.h
 isodate.o: isodate.h
 jacobi.o: Context.h Matrix.h blas.h
 jade.o: Context.h Matrix.h blas.h Timer.h
-kpgen.o: D3vector.h UnitCell.h
+kpgen.o: D3vector.h
 KpointCmd.o: UserInterface.h D3vector.h Sample.h AtomSet.h Context.h Atom.h
 KpointCmd.o: UnitCell.h ConstraintSet.h Wavefunction.h SharedFilePtr.h
 KpointCmd.o: Control.h
@@ -436,12 +437,12 @@ qb.o: SharedFilePtr.h Control.h Timer.h AngleCmd.h AtomCmd.h ComputeMLWFCmd.h
 qb.o: MLWFTransform.h BasisMapping.h ConstraintCmd.h DistanceCmd.h
 qb.o: FoldInWsCmd.h HelpCmd.h KpointCmd.h ListAtomsCmd.h ListSpeciesCmd.h
 qb.o: LoadCmd.h MoveCmd.h PrintCmd.h QuitCmd.h RandomizeWfCmd.h ResetVcmCmd.h
-qb.o: RunCmd.h SaveCmd.h SetCmd.h SpeciesCmd.h StatusCmd.h TorsionCmd.h
-qb.o: AtomsDyn.h Cell.h CellDyn.h SlaterDet.h Basis.h Matrix.h CellLock.h
-qb.o: CellMass.h ChargeMixCoeff.h ChargeMixRcut.h Debug.h Ecut.h Ecutprec.h
-qb.o: Ecuts.h Emass.h ExtStress.h FermiTemp.h Dt.h Nempty.h NetCharge.h
-qb.o: Nrowmax.h RefCell.h Stress.h Thermostat.h ThTemp.h ThTime.h ThWidth.h
-qb.o: WfDiag.h WfDyn.h Xc.h
+qb.o: RunCmd.h SaveCmd.h SetCmd.h SpeciesCmd.h StatusCmd.h StrainCmd.h
+qb.o: TorsionCmd.h AtomsDyn.h Cell.h CellDyn.h SlaterDet.h Basis.h Matrix.h
+qb.o: CellLock.h CellMass.h ChargeMixCoeff.h ChargeMixRcut.h Debug.h Ecut.h
+qb.o: Ecutprec.h Ecuts.h Emass.h ExtStress.h FermiTemp.h Dt.h Nempty.h
+qb.o: NetCharge.h Nrowmax.h RefCell.h Stress.h Thermostat.h ThTemp.h ThTime.h
+qb.o: ThWidth.h WfDiag.h WfDyn.h Xc.h
 qbox_xmlns.o: qbox_xmlns.h
 QuitCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
 QuitCmd.o: UnitCell.h ConstraintSet.h Wavefunction.h SharedFilePtr.h
@@ -529,6 +530,9 @@ spline.o: spline.h
 StatusCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
 StatusCmd.o: UnitCell.h ConstraintSet.h Wavefunction.h SharedFilePtr.h
 StatusCmd.o: Control.h
+StrainCmd.o: UserInterface.h Sample.h AtomSet.h Context.h Atom.h D3vector.h
+StrainCmd.o: UnitCell.h ConstraintSet.h Wavefunction.h SharedFilePtr.h
+StrainCmd.o: Control.h
 Stress.o: Sample.h AtomSet.h Context.h Atom.h D3vector.h UnitCell.h
 Stress.o: ConstraintSet.h Wavefunction.h SharedFilePtr.h Control.h
 StructuredDocumentHandler.o: StrX.h StructureHandler.h
@@ -622,10 +626,7 @@ XCPotential.o: BLYPFunctional.h Basis.h D3vector.h UnitCell.h
 XCPotential.o: FourierTransform.h blas.h
 XCPotential.o: ChargeDensity.h Timer.h Context.h LDAFunctional.h
 XCPotential.o: XCFunctional.h PBEFunctional.h BLYPFunctional.h
-xml2xy.o: StrX.h
 xmlextract.o: StrX.h
-xmlget.o: SampleHandler.h StructureHandler.h StructuredDocumentHandler.h
-xmlget.o: StrX.h
 XMLGFPreprocessor.o: Timer.h Context.h Base64Transcoder.h Matrix.h
 XMLGFPreprocessor.o: XMLGFPreprocessor.h
 XMLGFPreprocessor.o: Matrix.h

src/SDCellStepper.h

@@ -15,7 +15,7 @@
 // SDCellStepper.h:
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: SDCellStepper.h,v 1.5 2008-09-08 15:56:19 fgygi Exp $
+// $Id: SDCellStepper.h,v 1.6 2008-09-15 14:58:52 fgygi Exp $
 
 #ifndef SDCELLSTEPPER_H
 #define SDCELLSTEPPER_H
@@ -30,7 +30,7 @@ class SDCellStepper : public CellStepper
 
   SDCellStepper(Sample& s) : CellStepper(s) {}
 
-  void compute_new_cell(const std::valarray<double>& sigma_eks);
+  void compute_new_cell(const std::valarray<double>& sigma);
   void update_cell(void);
   double ekin(void) const { return 0.0; }
 };

src/StrainCmd.h

+////////////////////////////////////////////////////////////////////////////////
+//
+// 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:
+//
+////////////////////////////////////////////////////////////////////////////////
+// $Id: StrainCmd.h,v 1.1 2008-09-15 15:00:31 fgygi Exp $
+
+#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) {};
+
+  char *name(void) const { return "strain"; }
+  char *help_msg(void) const
+  {
+    return
+    "\n strain\n\n"
+    " syntax: strain [-atomsonly] [-inverse] uxx uxy uzz uxy uyz uxz \n\n"
+    "   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 =
+    " use: strain [-atomsonly] [-inverse] uxx uxy uzz uxy uyz uxz";
+    // 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

src/notes

 TODO: change all char* name() functions to const char* name ()
 (causes warning at compile time on tribe).
 --------------------------------------------------------------------------------
+rel1_45_0
+Added strain command for deformation of the cell.
+UnitCell.C: removed member compute_deda (moved calc to SDCellStepper). Added
+symmetric matrix - vector product function.
+SDCellStepper.C fixed calculation of cell correction from stress tensor.
+BOSampleStepper.C: reintroduced wf extrapolation when nempty>0 for MD.
+RunCmd.h: fixed help msg.
+--------------------------------------------------------------------------------
 rel1_44_5
 BOSampleStepper.C: Fix uninitialized mlwft ptr. Fixed delete statements.
 --------------------------------------------------------------------------------

src/qb.C

@@ -15,7 +15,7 @@
 // qb.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: qb.C,v 1.60 2008-09-08 15:56:20 fgygi Exp $
+// $Id: qb.C,v 1.61 2008-09-15 15:00:30 fgygi Exp $
 
 #include <iostream>
 #include <string>
@@ -62,6 +62,7 @@ using namespace std;
 #include "SetCmd.h"
 #include "SpeciesCmd.h"
 #include "StatusCmd.h"
+#include "StrainCmd.h"
 #include "TorsionCmd.h"
 
 #include "AtomsDyn.h"
@@ -246,6 +247,7 @@ int main(int argc, char **argv, char **envp)
   ui.addCmd(new SetCmd(s));
   ui.addCmd(new SpeciesCmd(s));
   ui.addCmd(new StatusCmd(s));
+  ui.addCmd(new StrainCmd(s));
   ui.addCmd(new TorsionCmd(s));
 
   ui.addVar(new AtomsDyn(s));

src/release.C

@@ -15,10 +15,10 @@
 // release.C
 //
 ////////////////////////////////////////////////////////////////////////////////
-// $Id: release.C,v 1.64 2008-09-09 04:25:39 fgygi Exp $
+// $Id: release.C,v 1.65 2008-09-15 15:00:30 fgygi Exp $
 
 #include "release.h"
 std::string release(void)
 {
-  return std::string("1.44.5");
+  return std::string("1.45.0");
 }