Added the complex absorbing potential in DVR calculations

Author: alexvakimov

src/dyn/wfcgrid2/Wfcgrid2.cpp

@@ -114,6 +114,7 @@ void Wfcgrid2::allocate(){
 
   Hdia = vector<CMATRIX>(Npts, CMATRIX(nstates, nstates));   ///<  diabatic Hamiltoninans for all the Npts points
   Hadi = vector<CMATRIX>(Npts, CMATRIX(nstates, nstates));   ///<  adiabatic Hamiltoninans for all the Npts points
+  Vcomplex = vector<CMATRIX>(Npts, CMATRIX(nstates, nstates));   ///<  complex absorbing potential
   U = vector<CMATRIX>(Npts, CMATRIX(nstates, nstates));      ///<  |adi> = |dia> * U : diabatic-to-adiabatic transformation for all the Npts points
   expH = vector<CMATRIX>(Npts, CMATRIX(nstates, nstates));   ///<  exponent of the diabatic Hamiltoninans for all the Npts points
   expK = vector<CMATRIX>(Npts, CMATRIX(nstates, nstates));   ///<  exponent of the kinetic energy propagator for all the Npts points

src/dyn/wfcgrid2/Wfcgrid2.h

@@ -134,6 +134,7 @@ class Wfcgrid2{
   /// Hamiltonian and Propagator
   vector<CMATRIX> Hdia;      ///<  diabatic Hamiltoninans for all the Npts points
   vector<CMATRIX> Hadi;      ///<  adiabatic Hamiltoninans for all the Npts points
+  vector<CMATRIX> Vcomplex;  ///<  complex absorbing potential that will be added to the real-space propagator
   vector< vector<CMATRIX> > NAC1;  ///<  1-st order NACS: NAC1[ipt][alpha].get(i,j) = <psi_i| nabla_alpha | psi_j> 
   vector< vector<CMATRIX> > NAC2;  ///<  2-nd order NACS: NAC2[ipt][alpha].get(i,j) = <psi_i| nabla_alpha^2 | psi_j> 
   vector<CMATRIX> U;         ///<  |adi> = |dia> * U : diabatic-to-adiabatic transformation for all the Npts points
@@ -231,6 +232,7 @@ class Wfcgrid2{
   ///< Precompute Hamiltonian on the grid
   void update_Hamiltonian(bp::object py_funct, bp::object params, int rep);
 
+
   ///< Convert diabatic and adiabatic wfc into one another
   void update_adiabatic();
   void update_diabatic();

src/dyn/wfcgrid2/Wfcgrid2_SOFT.cpp

@@ -36,40 +36,59 @@ using namespace libmeigen;
 
 void Wfcgrid2::update_propagator_H(double dt){
 /**
-  \brief Update the Hamiltonian for nd-D grid
+  \brief Update the exp(-i*dt*Hamiltonian) for nd-D grid
 
   Should have called ```update_Hamiltonian``` prior to this
 
+  H * U = U * E  =>  H = U * E * U^+
+
+  exp(-i*dt*H) = U * exp(-i*dt*E) * U^+
+
+  exp(-i*dt*E) = exp(-i * dt * [Re(E)+i*Im(E)]) = exp(-i*dt*Re(E)) * exp(dt*Im(E))
+                                                  ----- ex -------   ---- scl ----
+
+  ex = cos(dt*Re(E)) - i * sin(dt*Re(E))
 */
 
   int nst, nst1;
 
   CMATRIX S(nstates, nstates); S.Init_Unit_Matrix(1.0);
-  CMATRIX si(nstates, nstates);  // sin(-dt*E), where E is adiabatic Ham.
-  CMATRIX cs(nstates, nstates);  // cos(-dt*E), where E is adiabatic Ham.
+//  CMATRIX si(nstates, nstates);  // sin(dt*E), where E is adiabatic Ham.
+//  CMATRIX cs(nstates, nstates);  // cos(dt*E), where E is adiabatic Ham.
+//  CMATRIX ex(nstates, nstates);  // exp()
 
   complex<double> one(1.0, 0.0);
   complex<double> eye(0.0, 1.0);
+  //complex<double> ex(0.0, 0.0);
+  double scl = 1.0;
+  double si, cs;
 
   // For all grid points
   for(int npt1=0; npt1<Npts; npt1++){
 
     // Transformation to adiabatic basis
     // Hdia * U = S * U * Hadi
     solve_eigen(Hdia[npt1], S, Hadi[npt1], U[npt1], 0);  
-
-    cs = 0.0;   si = 0.0;
+    
+//    cs = 0.0;   si = 0.0;
     for(int nst=0;nst<nstates;nst++){  
-      cs.set(nst, nst, std::cos(-dt*Hadi[npt1].get(nst, nst)) * one );
-      si.set(nst, nst, std::sin(-dt*Hadi[npt1].get(nst, nst)) * one );
+      //cs.set(nst, nst, std::cos(dt*Hadi[npt1].get(nst, nst)) * one );
+      //si.set(nst, nst, std::sin(dt*Hadi[npt1].get(nst, nst)) * one );
+      cs = std::cos( dt*( Hadi[npt1].get(nst, nst).real() + Vcomplex[npt1].get(nst, nst).real() ) );
+      si = std::sin( dt*( Hadi[npt1].get(nst, nst).real() + Vcomplex[npt1].get(nst, nst).real() ) );
+      scl = std::exp(dt*( Hadi[npt1].get(nst, nst).imag() + Vcomplex[npt1].get(nst, nst).imag() ) );
+       
+      complex<double> ex(cs, -si);
+      expH[npt1].set(nst, nst,  scl * ex);
     }
 
     // Transform cs and si according to matrix U:
-    cs = U[npt1] * cs * U[npt1].H(); 
-    si = U[npt1] * si * U[npt1].H(); 
+    //cs = U[npt1] * cs * U[npt1].H(); 
+    //si = U[npt1] * si * U[npt1].H(); 
+    expH[npt1] = U[npt1] * expH[npt1] * U[npt1].H();
 
     // Finally construct complex exp(-i*dt*H) matrix from the cs and si matrices
-    expH[npt1] = cs + eye*si;
+    //expH[npt1] = cs - eye*si;
 
 
   }// for npt1  
@@ -78,6 +97,7 @@ void Wfcgrid2::update_propagator_H(double dt){
 
 
 
+
 void Wfcgrid2::update_propagator_K(double dt, vector<double>& mass){
 /**
   \brief Update reciprocal-space propagators for nd-D grid

src/dyn/wfcgrid2/Wfcgrid2_updates.cpp

@@ -1,8 +1,8 @@
 /*********************************************************************************
-* Copyright (C) 2019 Alexey V. Akimov
+* Copyright (C) 2019-2022 Alexey V. Akimov
 *
 * This file is distributed under the terms of the GNU General Public License
-* as published by the Free Software Foundation, either version 2 of
+* as published by the Free Software Foundation, either version 3 of
 * the License, or (at your option) any later version.
 * See the file LICENSE in the root directory of this distribution
 * or <http://www.gnu.org/licenses/>.
@@ -58,6 +58,11 @@ void Wfcgrid2::update_Hamiltonian(bp::object py_funct, bp::object params, int re
     bp::object obj = py_funct(bp::object(q), params);  
 
 
+    // Try extract the adiabatic Hamiltonian
+    has_attr = (int)hasattr(obj,"v_complex");
+    if(has_attr){   Vcomplex[npt1] = extract<CMATRIX>(obj.attr("v_complex"));    } 
+
+
     if(rep==0){
 
       // Try to extract the diabatic Hamiltonian
@@ -90,6 +95,8 @@ void Wfcgrid2::update_Hamiltonian(bp::object py_funct, bp::object params, int re
 
 
 
+
+
 void Wfcgrid2::update_adiabatic(){
 /**
   Update the adiabatic wfc from the diabatic: DIA -> ADI

src/dyn/wfcgrid2/libwfcgrid2.cpp

@@ -151,6 +151,7 @@ void export_Wfcgrid2_objects(){
 
       .def_readwrite("Hdia", &Wfcgrid2::Hdia)
       .def_readwrite("Hadi", &Wfcgrid2::Hadi)
+      .def_readwrite("Vcomplex", &Wfcgrid2::Vcomplex)
       .def_readwrite("NAC1", &Wfcgrid2::NAC1)
       .def_readwrite("NAC2", &Wfcgrid2::NAC2)
       .def_readwrite("U",    &Wfcgrid2::U)