Cleaned up a bit some of the F-tracking procedures, made them more efficient, fixed some exports, added some experimental versions of the F-tracking, but not what's the current manuscript uses - something to consider in the future

Author: alexvakimov

src/dyn/dyn_ham.cpp

@@ -71,6 +71,7 @@ void update_Hamiltonian_variables(dyn_control_params& prms, dyn_variables& dyn_v
 */ 
 
   int nadi = ham.nadi;
+  int ndof = dyn_var.ndof;
   int ntraj = ham.children.size();
 
   MATRIX& q = *dyn_var.q;
@@ -147,6 +148,44 @@ void update_Hamiltonian_variables(dyn_control_params& prms, dyn_variables& dyn_v
 
   if(update_type==0 || update_type==1){
 
+
+    // Derivative NAC correction option:
+    if(prms.do_nac_phase_correction==2){  // Experimental option to fix the phase of NACVs
+
+      vector<int> traj_id(1,0);
+      for(int traj=0; traj<ntraj; traj++){
+        traj_id[0] = traj;
+        CMATRIX Eadi(ham.children[traj]->get_ham_adi());
+        MATRIX e_curr(ham.children[traj]->get_ham_adi().real());
+        MATRIX e_prev(ham_prev.children[traj]->get_ham_adi().real());
+        CMATRIX f_curr = ham.children[traj]->all_forces_adi(traj_id);
+        CMATRIX f_prev = ham_prev.children[traj]->all_forces_adi(traj_id);
+        //ham->get_hvib_adi().show_matrix();
+        MATRIX pp(dyn_var.ndof, 1); 
+        double dt = prms.dt;
+        pp = p.col(traj);
+        int act_state = dyn_var.act_states[traj];
+
+        //vector<MATRIX> T_new(compute_F_cost_matrix_dof_resolved(f_curr, f_prev, e_curr, e_prev, pp, iM, dt, act_state));
+        MATRIX T_new(compute_F_cost_matrix(f_curr, f_prev, e_curr, e_prev, pp, iM, dt, act_state).real() );
+
+        for(int k=0; k<ndof;k++){
+          for(int i=0; i<nadi;i++){
+            for(int j=i+1; j<nadi; j++){
+              //double sgn = T_new[k].get(i,j);
+              double sgn = T_new.get(i,j);
+              if(sgn > 0.5){ sgn = -1; } // change sign
+              else{ sgn = 1.0; }
+              
+              ham.children[traj]->dc1_adi[k]->scale(i,j, sgn);
+              ham.children[traj]->dc1_adi[k]->scale(j,i, sgn);
+            }// for j
+          }// for i
+
+        }// for k
+      }// for traj
+    }// if correction
+
     //========================== Couplings ===============================
     // Don't update NACs - they may have been read in step 1
     if(prms.nac_update_method==0){ ;;  }    

src/dyn/dyn_projectors.cpp

@@ -1115,60 +1115,111 @@ CMATRIX compute_F_cost_matrix(CMATRIX& F_curr,  CMATRIX& F_prev, MATRIX& e_curr,
 
   MATRIX vel(ndof, 1);  vel.dot_product(iM, p);  
   MATRIX tmp(ndof, 1);
-  tmp = fi_curr.T(); ai.dot_product(iM, tmp);
-  tmp = fj_curr.T(); aj.dot_product(iM, tmp); 
 
   double val = 0.0;
   for(i=0; i<nst; i++){
-    double sum = 1.0; 
 
-    for(j=0;j<nst; j++){
+    fi_curr = F_curr.row(i).real();
+    fi_prev = F_prev.row(i).real();
+    tmp = fi_curr.T(); ai.dot_product(iM, tmp);
+    dq_i = (vel * dt + 0.5 * ai * dt * dt);
+    double dE_i = -(fi_prev * dq_i).get(0,0);
 
-      fi_curr = F_curr.row(i).real();
-      fi_prev = F_prev.row(i).real();
+    for(j=0;j<nst; j++){
 
       fj_curr = F_curr.row(j).real();
       fj_prev = F_prev.row(j).real();
-
-      dq_i = (vel * dt + 0.5 * ai * dt * dt); 
-      double dE = e_prev.get(j,j) - e_prev.get(i,i);
-      double dE_i = -(fi_prev * dq_i).get(0,0);
-
+      tmp = fj_curr.T(); aj.dot_product(iM, tmp);
       dq_j = (vel * dt + 0.5 * aj * dt * dt);
       double dE_j = -(fj_prev * dq_j).get(0,0);
+
+      double dE = e_prev.get(j,j) - e_prev.get(i,i);
       double dE_new = dE + dE_j - dE_i;
 
       val = 0.5*(1.0 - SIGN( dE ) * SIGN( dE_new )) ;
-/*
-      if(i==j){ val = 0.0; }
-      //else{ val = 1.0 - SIGN( dE ) * SIGN( dE_new );  sum -= val; } 
-      else{ 
-         double val1;
-         val = SIGN( dE ) * SIGN( dE_new );  
-         if(val<0.0){ val1 = 1.0;  }
-         else{ val1 = 0.0; }
-         val = val1;
-         //sum -= val;
-      }
-*/
 
       res.set(i, j, complex<double>( val, 0.0) );
     }// for j
-
-    //if(sum<=0.0){ sum = 0.0; }
-    //res.set(i, i, complex<double>(sum, 0.0) );
   }// for i
 
   return res;
 }
 
 
+
+vector<MATRIX> compute_F_cost_matrix_dof_resolved(CMATRIX& F_curr,  CMATRIX& F_prev, MATRIX& e_curr, MATRIX& e_prev, MATRIX& _p, MATRIX& iM, double dt, int act_state){
+/**  F_curr, F_prev - CMATRIX(nadi, ndof)
+     p - MATRIX(ndof, 1) - current momentum
+     iM - MATRIX(ndof, 1) - inverse matrices
+     dt - time step
+     act_state - the current active state
+*/
+  int i,j;
+  int ndof = F_curr.n_cols;
+  int nst = F_curr.n_rows;
+  vector<MATRIX> res(ndof, MATRIX(nst, nst) );
+  //MATRIX fi_curr(1, ndof);
+  //MATRIX fj_curr(1, ndof);
+  //MATRIX fi_prev(1, ndof);
+  //MATRIX fj_prev(1, ndof);
+  //MATRIX ai(ndof, 1);
+  //MATRIX aj(ndof, 1);
+  //MATRIX dq_i(ndof,1);
+  //MATRIX dq_j(ndof,1);
+  MATRIX p(_p);
+
+  // correction, because the current p is already a half-step propagated using old forces
+  //p -= dt * F_prev.row(act_state).real().T();
+
+  MATRIX vel(ndof, 1);  vel.dot_product(iM, p);
+  MATRIX tmp(ndof, 1);
+
+  double val = 0.0;
+  for(int idof=0; idof<ndof; idof++){
+
+    double vel = p.get(idof, 0) * iM.get(idof, 0);
+
+    for(i=0; i<nst; i++){
+
+      double fi_curr = F_curr.get(i,idof).real();
+      double fi_prev = F_prev.get(i,idof).real();
+      double ai = fi_curr * iM.get(idof,0);
+      double dq_i = vel * dt + 0.5 * ai * dt * dt;
+      double dE_i = -fi_prev * dq_i;
+
+      for(j=0;j<nst; j++){
+
+        double fj_curr = F_curr.get(j,idof).real();
+        double fj_prev = F_prev.get(j,idof).real();
+        double aj = fj_curr * iM.get(idof,0);
+        double dq_j = vel * dt + 0.5 * ai * dt * dt;
+        double dE_j = -fj_prev * dq_j;
+
+        double dE = e_prev.get(j,j) - e_prev.get(i,i);
+        double dE_new = dE + dE_j - dE_i;
+
+        val =  SIGN( dE ) * SIGN( dE_new ) ;
+        res[idof].set(i, j, val);
+
+      }// for j
+    }// for i
+  }// for idof
+
+  return res;
+}
+
+
+
+
 CMATRIX compute_F_cost_matrix2(CMATRIX& F_curr,  CMATRIX& F_prev, MATRIX& e_curr, MATRIX& e_prev, MATRIX& _p, MATRIX& iM, double dt, int act_state){
 /**  F_curr, F_prev - CMATRIX(nadi, ndof)
      p - MATRIX(ndof, 1) - current momentum
      iM - MATRIX(ndof, 1) - inverse matrices
      dt - time step
      act_state - the current active state
+
+   This is an experimental function - merely a placeholder - not used anywhere yet!
+
 */
   int i,j;
   int ndof = F_curr.n_cols;

src/dyn/dyn_projectors.h

@@ -64,6 +64,9 @@ CMATRIX raw_to_dynconsyst(CMATRIX& amplitudes, vector<CMATRIX>& projectors);
 CMATRIX dynconsyst_to_raw(CMATRIX& amplitudes, vector<CMATRIX>& projectors);
 
 CMATRIX compute_F_cost_matrix(CMATRIX& F_curr,  CMATRIX& F_prev, MATRIX& e_curr, MATRIX& e_prev, MATRIX& p, MATRIX& iM, double dt, int act_state);
+vector<MATRIX> compute_F_cost_matrix_dof_resolved(CMATRIX& F_curr,  CMATRIX& F_prev, MATRIX& e_curr, MATRIX& e_prev, 
+MATRIX& _p, MATRIX& iM, double dt, int act_state);
+
 
 }// namespace libdyn
 }// liblibra

src/dyn/libdyn.cpp

@@ -802,6 +802,15 @@ void export_dyn_projectors_objects(){
   def("get_stochastic_reordering3", expt_get_stochastic_reordering3_v2);
 
 
+  CMATRIX (*expt_compute_F_cost_matrix_v1)
+  (CMATRIX& F_curr,  CMATRIX& F_prev, MATRIX& e_curr, MATRIX& e_prev, 
+  MATRIX& p, MATRIX& iM, double dt, int act_state) = &compute_F_cost_matrix;
+  def("compute_F_cost_matrix", expt_compute_F_cost_matrix_v1);
+
+  vector<MATRIX> (*expt_compute_F_cost_matrix_dof_resolved_v1)
+  (CMATRIX& F_curr,  CMATRIX& F_prev, MATRIX& e_curr, MATRIX& e_prev,
+  MATRIX& _p, MATRIX& iM, double dt, int act_state) = &compute_F_cost_matrix_dof_resolved;
+  def("compute_F_cost_matrix_dof_resolved", expt_compute_F_cost_matrix_dof_resolved_v1);
 
 }