Added kpoint move functionality

src/KpointCmd.h

@@ -39,6 +39,7 @@ class KpointCmd : public Cmd
     " syntax:\n\n"
     "   kpoint add kx ky kz weight\n"
     "   kpoint delete kx ky kz \n"
+    "   kpoint move kx ky kz new_kx new_ky new_kz\n"
     "   kpoint list\n\n";
   }
 
@@ -79,6 +80,22 @@ class KpointCmd : public Cmd
       double kz = atof(argv[4]);
       s->wf.del_kpoint(D3vector(kx,ky,kz));
     }
+    else if ( subcmd == "move" )
+    {
+      if ( argc != 8 )
+      {
+        if ( onpe0 )
+          cout << help_msg();
+        return 1;
+      }
+      double kx = atof(argv[2]);
+      double ky = atof(argv[3]);
+      double kz = atof(argv[4]);
+      double newkx = atof(argv[5]);
+      double newky = atof(argv[6]);
+      double newkz = atof(argv[7]);
+      s->wf.move_kpoint(D3vector(kx,ky,kz),D3vector(newkx,newky,newkz));
+    }
     else if ( subcmd == "list" )
     {
       if ( argc != 2 )
@@ -89,7 +106,7 @@ class KpointCmd : public Cmd
       }
       if ( onpe0 )
       {
-        cout << " <-- kpoint list: reciprocal lattice coordinates" << endl;
+        cout << " <!-- kpoint list: reciprocal lattice coordinates" << endl;
         for ( int ikp = 0; ikp < s->wf.nkp(); ikp++ )
         {
           cout << " "

src/Wavefunction.C

@@ -18,6 +18,7 @@
 
 #include "Wavefunction.h"
 #include "SlaterDet.h"
+#include "FourierTransform.h"
 #include "jacobi.h"
 #include "SharedFilePtr.h"
 #include <vector>
@@ -357,7 +358,7 @@ void Wavefunction::add_kpoint(D3vector kpoint, double weight)
     {
       if ( ctxt_.onpe0() )
         cout << " Wavefunction::add_kpoint: kpoint already defined"
-           << endl;
+             << endl;
       return;
     }
   }
@@ -426,6 +427,108 @@ void Wavefunction::del_kpoint(D3vector kpoint)
 }
 
 ////////////////////////////////////////////////////////////////////////////////
+void Wavefunction::move_kpoint(D3vector kpoint, D3vector newkpoint)
+{
+  bool found = false;
+  int ikp = 0;
+  while ( !found && ikp < kpoint_.size() )
+  {
+    if ( length(kpoint_[ikp] - kpoint) < 1.e-6 )
+    {
+      found = true;
+    }
+    else
+    {
+      ikp++;
+    }
+  }
+  if ( !found )
+  {
+    if ( ctxt_.onpe0() )
+      cout << " Wavefunction::move_kpoint: no such kpoint"
+         << endl;
+    return;
+  }
+  // check if newkpoint coincides with existing kpoint
+  for ( int i = 0; i < kpoint_.size(); i++ )
+  {
+    if ( length(newkpoint - kpoint_[i]) < 1.e-6 )
+    {
+      if ( ctxt_.onpe0() )
+        cout << " Wavefunction::move_kpoint: kpoint already defined "
+             << "at newkpoint position"
+             << endl;
+      return;
+    }
+  }
+
+  // copy wavefunctions from old SlaterDet at kpoint to new SlaterDet
+  // at newkpoint
+  for ( int ispin = 0; ispin < nspin_; ispin++ )
+  {
+    // create new SlaterDet at newkpoint
+    SlaterDet *sd = sd_[ispin][ikp];
+    SlaterDet *sdn = new SlaterDet(*sdcontext_,newkpoint);
+    sdn->resize(cell_,refcell_,ecut_,nst_[ispin]);
+    sdn->init();
+    // copy wave functions from old to new SlaterDet
+    const Basis& basis = sd_[ispin][ikp]->basis();
+    const Basis& newbasis = sdn->basis();
+    // transform all states to real space and interpolate
+    int np0 = max(basis.np(0),newbasis.np(0));
+    int np1 = max(basis.np(1),newbasis.np(1));
+    int np2 = max(basis.np(2),newbasis.np(2));
+    FourierTransform ft1(basis,np0,np1,np2);
+    FourierTransform ft2(newbasis,np0,np1,np2);
+    // allocate real-space grid
+    valarray<complex<double> > tmpr(ft1.np012loc());
+    for ( int n = 0; n < sd->nstloc(); n++ )
+    {
+      for ( int i = 0; i < tmpr.size(); i++ )
+        tmpr[i] = 0.0;
+      ComplexMatrix& c = sd->c();
+      ComplexMatrix& cn = sdn->c();
+      ft1.backward(c.valptr(n*c.mloc()),&tmpr[0]);
+      // if the new kpoint is Gamma, remove the phase of the wavefunction
+      // using the phase of the first element of the array
+      if ( newbasis.real() )
+      {
+        double arg0 = arg(tmpr[0]);
+        // broadcast argument found on task 0
+        MPI_Bcast(&arg0,1,MPI_DOUBLE,0,basis.comm());
+        complex<double> z = polar(1.0,-arg0);
+        for ( int i = 0; i < tmpr.size(); i++ )
+          tmpr[i] *= z;
+      }
+      ft2.forward(&tmpr[0], cn.valptr(n*cn.mloc()));
+    }
+
+    // delete old SlaterDet
+    delete sd_[ispin][ikp];
+    // reassign pointer
+    sd_[ispin][ikp] = sdn;
+  }
+  kpoint_[ikp] = newkpoint;
+
+  if ( nspin_ == 1 )
+  {
+    sd_[0][ikp]->update_occ(nel_,nspin_);
+  }
+  else if ( nspin_ == 2 )
+  {
+    const int nocc_up = (nel_+1)/2+deltaspin_;
+    const int nocc_dn = nel_/2 - deltaspin_;
+    sd_[0][ikp]->update_occ(nocc_up,nspin_);
+    sd_[1][ikp]->update_occ(nocc_dn,nspin_);
+  }
+  else
+  {
+    // incorrect value of nspin_
+    assert(false);
+  }
+}
+
+////////////////////////////////////////////////////////////////////////////////
 void Wavefunction::randomize(double amplitude)
 {
   for ( int ispin = 0; ispin < nspin_; ispin++ )

src/Wavefunction.h

@@ -101,6 +101,7 @@ class Wavefunction
   void set_nrowmax(int n);
   void add_kpoint(D3vector kpoint, double weight);
   void del_kpoint(D3vector kpoint);
+  void move_kpoint(D3vector kpoint, D3vector new_kpoint);
 
   void randomize(double amplitude);