Merge branch 'devel'

Author: alexvakimov

CMakeLists.txt

@@ -96,9 +96,12 @@ SET(EIGEN3_INCLUDE_DIRS "${EIGEN3_INCLUDE_DIR}")
 
 # Libint2 2.7.0 CONFIG REQUIRED COMPONENTS shared gss e5 g5
 #FIND_PACKAGE(Libint2 MODULE 2.7.0)
-FIND_PACKAGE(Libint2 CONFIG 2.7.0 REQUIRED COMPONENTS shared gss e5 g5)
+#find_package(Libint2)
+FIND_PACKAGE(Libint2 CONFIG 2.7.1 REQUIRED COMPONENTS shared gss impure_sh onebody_d0_l6 g12_d0_l4 g12_d1_l4 eri_c4_d0_l5 eri_c4_d1_l4) #gss e5 g5)
 IF(Libint2_FOUND)
-  INCLUDE_DIRECTORIES("${Libint2_INCLUDE_DIRS}")   
+  #INCLUDE_DIRECTORIES("${Libint2_INCLUDE_DIRS}")   
+  INCLUDE_DIRECTORIES(Libint2::int2 Libint2::cxx Libint2::int2-cxx Libint2::impure_sh Libint2::gss Libint2::onebody_d0_l6 Libint::g12_d0_l4 Libint2::g12_d1_l4 Libint2::shared Libint2::cxx_ho Libint2::c  Libint2::eri_c4_d0_l5  Libint2::eri_c4_d1_l4)   
+  MESSAGE("${LIBINT2_INCLUDE_DIR}")
   MESSAGE("Found Libint2 include directory: ")
   MESSAGE("${Libint2_INCLUDE_DIRS}")
   MESSAGE("Found Libint2_LIBRARIES: ")
@@ -136,7 +139,7 @@ ENDIF()
 #
 # Set the libraries
 # 
-SET( ext_libs ${Boost_LIBRARIES} ${PYTHON_LIBRARIES} ) 
+SET( ext_libs ${Boost_LIBRARIES} ${PYTHON_LIBRARIES}) 
 
 
 

README.md

@@ -26,6 +26,20 @@ for all users with the intent:
    knowledge and skills with others;
 
 
+## Installation Videotutorials (as of 5/16/2022)
+
+* [Installing WSL2](https://ub.hosted.panopto.com/Panopto/Pages/Embed.aspx?id=02184b70-7745-4eb4-a776-ae92014c652a&autoplay=false&offerviewer=true&showtitle=true&showbrand=true&captions=false&interactivity=all)
+
+* [Installing WSL2: After reboot](https://ub.hosted.panopto.com/Panopto/Pages/Embed.aspx?id=972aef79-e235-4a90-9ce1-ae92014d34db&autoplay=false&offerviewer=true&showtitle=true&showbrand=true&captions=false&interactivity=all)
+
+* [Installing Ubuntu on Windows 11](https://ub.hosted.panopto.com/Panopto/Pages/Embed.aspx?id=31a63536-f333-4242-9b56-ae92015ece64&autoplay=false&offerviewer=true&showtitle=true&showbrand=true&captions=false&interactivity=all)
+
+* [Creating Conda environment for Libra](https://ub.hosted.panopto.com/Panopto/Pages/Embed.aspx?id=d6ada23e-7e16-4b7a-b290-ae920188627c&autoplay=false&offerviewer=true&showtitle=true&showbrand=true&captions=false&interactivity=all)
+
+* [Installing Libra](https://ub.hosted.panopto.com/Panopto/Pages/Embed.aspx?id=7f8dd8c4-9f58-4ca0-a8cb-ae930166b7ec&autoplay=false&offerviewer=true&showtitle=true&showbrand=true&captions=false&interactivity=all)
+
+
+
 ## Installation (as of after 5/14/2021)
 
 ### 1. Install miniconda (for Python 3.8) and activate Conda
@@ -131,7 +145,7 @@ We need to downgrade Python version here to 3.6 to enable Psi4 installation
 You can install Jupyter notebook using the following command. This will be useful and you can set up and access the Jupyter notebook 
 remotely from a cluster and load the tutorials
 
-    conda install -c conda-forge jupyterlab
+    conda install -c anaconda jupyter
 
 
 Used in some of the tutorials

src/dyn/Dynamics.cpp

@@ -248,7 +248,7 @@ CMATRIX transform_amplitudes(int rep_in, int rep_out, CMATRIX& C, nHamiltonian&
 
 
 //vector<CMATRIX> compute_St(nHamiltonian& ham, CMATRIX** Uprev){
-vector<CMATRIX> compute_St(nHamiltonian& ham, vector<CMATRIX>& Uprev){
+vector<CMATRIX> compute_St(nHamiltonian& ham, vector<CMATRIX>& Uprev, int isNBRA){
 /**
   This function computes the time-overlap matrices for all trajectories
 
@@ -258,16 +258,26 @@ vector<CMATRIX> compute_St(nHamiltonian& ham, vector<CMATRIX>& Uprev){
   int ntraj = ham.children.size();
 
   vector<CMATRIX> St(ntraj, CMATRIX(nst, nst));
-
+  if(isNBRA==1){
+    St[0] = Uprev[0].H() * ham.children[0]->get_basis_transform();
+  }
+  else{
   for(int traj=0; traj<ntraj; traj++){
     St[traj] = Uprev[traj].H() * ham.children[traj]->get_basis_transform();
   }
-
+  }
   return St;
 
 }
 
-vector<CMATRIX> compute_St(nHamiltonian& ham){
+vector<CMATRIX> compute_St(nHamiltonian& ham, vector<CMATRIX>& Uprev){
+
+  int is_nbra = 0;
+  return compute_St(ham, Uprev, is_nbra);
+
+}
+
+vector<CMATRIX> compute_St(nHamiltonian& ham, int isNBRA){
 /**
   This function computes the time-overlap matrices for all trajectories
 
@@ -277,16 +287,23 @@ vector<CMATRIX> compute_St(nHamiltonian& ham){
   int ntraj = ham.children.size();
 
   vector<CMATRIX> St(ntraj, CMATRIX(nst, nst));
-
+  if(isNBRA==1){
+    St[0] = ham.children[0]->get_time_overlap_adi();
+  }
+  else{
   for(int traj=0; traj<ntraj; traj++){
     St[traj] = ham.children[traj]->get_time_overlap_adi();
   }
-
+  }
   return St;
 
 }
 
+vector<CMATRIX> compute_St(nHamiltonian& ham){
+  int is_nbra = 1;
 
+  return compute_St(ham, is_nbra);
+}
 
 
 void apply_afssh(dyn_variables& dyn_var, CMATRIX& C, vector<int>& act_states, MATRIX& invM,
@@ -487,20 +504,20 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
 
 
 
-
   dyn_control_params prms;
   prms.set_parameters(dyn_params);
 
   int i,j;
   int cdof;
   int ndof = q.n_rows;
   int ntraj = q.n_cols;
+  int ntraj1;
   int nst = C.n_rows;    
   int traj, dof, idof;
   int n_therm_dofs;
   int num_el = prms.num_electronic_substeps;
   double dt_el = prms.dt / num_el;
-
+  //cout << "Flag 1 isNBRA " << prms.isNBRA << endl;
 
   //dyn_variables dyn_var(nst, nst, ndof, ntraj);
 
@@ -510,11 +527,22 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
   vector<int> project_out_states(ntraj); // for DISH
 
   vector<CMATRIX> Uprev;
-  vector<CMATRIX> St(ntraj, CMATRIX(nst, nst));
-  vector<CMATRIX> Eadi(ntraj, CMATRIX(nst, nst));  
-  vector<MATRIX> decoherence_rates(ntraj, MATRIX(nst, nst)); 
-  vector<double> Ekin(ntraj, 0.0);  
-  vector<MATRIX> prev_ham_dia(ntraj, MATRIX(nst, nst));
+  // Defining ntraj1 as a reference for making these matrices
+  // ntraj is defined as q.n_cols as since it would be large in NBRA
+  // we can define another variable like ntraj1 and build the matrices based on that.
+  // We can make some changes where q is generated but this seems to be a bit easier
+  if(prms.isNBRA==1){
+  ntraj1 = 1;
+  }
+  else{
+  ntraj1 = ntraj;
+  }
+  // Defining matrices based on ntraj1
+  vector<CMATRIX> St(ntraj1, CMATRIX(nst, nst));
+  vector<CMATRIX> Eadi(ntraj1, CMATRIX(nst, nst));  
+  vector<MATRIX> decoherence_rates(ntraj1, MATRIX(nst, nst)); 
+  vector<double> Ekin(ntraj1, 0.0);  
+  vector<MATRIX> prev_ham_dia(ntraj1, MATRIX(nst, nst));
   MATRIX gamma(ndof, ntraj);
   MATRIX p_traj(ndof, 1);
   vector<int> t1(ndof, 0); for(dof=0;dof<ndof;dof++){  t1[dof] = dof; }
@@ -535,13 +563,14 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
   }
 
 
-
-  if(prms.tsh_method == 3){ // DISH
-    for(traj=0; traj<ntraj; traj++){
+ if(prms.tsh_method == 3){ // DISH
+  //    prev_ham_dia[0] = ham.children[0]->get_ham_dia().real();
+  //}
+  //else{
+    for(traj=0; traj<ntraj1; traj++){
       prev_ham_dia[traj] = ham.children[traj]->get_ham_dia().real();  
     }
   }
-
   //============ Update the Hamiltonian object =============
   // In case, we may need phase correction & state reordering
   // prepare the temporary files
@@ -562,13 +591,12 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
 
   //============== Electronic propagation ===================
   // Evolve electronic DOFs for all trajectories
+  // Adding the prms.isNBRA to the propagate electronic
   for(i=0; i<num_el; i++){
-    propagate_electronic(0.5*dt_el, C, projectors, ham.children, prms.rep_tdse);   
+    propagate_electronic(0.5*dt_el, C, projectors, ham.children, prms.rep_tdse, prms.isNBRA);   
   }
 
-
   //============== Nuclear propagation ===================
-
   // NVT dynamics
   if(prms.ensemble==1){  
     for(idof=0; idof<n_therm_dofs; idof++){
@@ -578,9 +606,7 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
       }// traj
     }// idof 
   }
- 
   p = p + aux_get_forces(prms, C, projectors, act_states, ham) * 0.5 * prms.dt;
-
   // Kinetic constraint
   for(cdof = 0; cdof < prms.constrained_dofs.size(); cdof++){   
     p.scale(prms.constrained_dofs[cdof], -1, 0.0); 
@@ -591,7 +617,6 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
   if(prms.entanglement_opt==22){
     gamma = ETHD3_friction(q, p, invM, prms.ETHD3_alpha, prms.ETHD3_beta);
   }
-
   // Update coordinates of nuclei for all trajectories
   for(traj=0; traj<ntraj; traj++){
     for(dof=0; dof<ndof; dof++){  
@@ -615,9 +640,10 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
 
     if(prms.state_tracking_algo > 0 || prms.do_phase_correction){
 
-      if(prms.time_overlap_method==0){    St = compute_St(ham, Uprev);    }
-      else if(prms.time_overlap_method==1){    St = compute_St(ham);      }
-
+      if(prms.time_overlap_method==0){    
+        St = compute_St(ham, Uprev, prms.isNBRA);
+      }
+      else if(prms.time_overlap_method==1){    St = compute_St(ham, prms.isNBRA);      }
       Eadi = get_Eadi(ham);           // these are raw properties
       update_projectors(prms, projectors, Eadi, St, rnd);
 
@@ -645,7 +671,6 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
     p.scale(prms.constrained_dofs[cdof], -1, 0.0); 
   }
 
-
   // NVT dynamics
   if(prms.ensemble==1){  
     for(idof=0; idof<n_therm_dofs; idof++){
@@ -656,13 +681,13 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
     }// idof 
   }
 
-
   //============== Electronic propagation ===================
   // Evolve electronic DOFs for all trajectories
   update_Hamiltonian_p(prms, ham, p, invM);
   for(i=0; i<num_el; i++){
-    propagate_electronic(0.5*dt_el, C, projectors, ham.children, prms.rep_tdse);   
+    propagate_electronic(0.5*dt_el, C, projectors, ham.children, prms.rep_tdse, prms.isNBRA);
   }
+  CMATRIX Hvib(ham.children[0]->nadi, ham.children[0]->nadi);  Hvib = ham.children[0]->get_hvib_adi();
 
 
   //============== Begin the TSH part ===================
@@ -684,27 +709,25 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
   }
 
 
-
   //================= Update decoherence rates & times ================
   //MATRIX decoh_rates(*prms.decoh_rates);
 //    exit(0);
-  
   if(prms.decoherence_times_type==-1){
-    for(traj=0; traj<ntraj; traj++){   decoherence_rates[traj] = 0.0;   }
+    for(traj=0; traj<ntraj1; traj++){   decoherence_rates[traj] = 0.0;   }
   }
 
   /// mSDM
   /// Just use the plain times given from the input, usually the
   /// mSDM formalism
   else if(prms.decoherence_times_type==0){
-    for(traj=0; traj<ntraj; traj++){   decoherence_rates[traj] = *prms.decoherence_rates;   }
+    for(traj=0; traj<ntraj1; traj++){   decoherence_rates[traj] = *prms.decoherence_rates;   }
   }
 
   /// Compute the dephasing rates according the original energy-based formalism
   else if(prms.decoherence_times_type==1){
     Eadi = get_Eadi(ham); 
     Ekin = compute_kinetic_energies(p, invM);
-    decoherence_rates = edc_rates(Eadi, Ekin, prms.decoherence_C_param, prms.decoherence_eps_param);       
+    decoherence_rates = edc_rates(Eadi, Ekin, prms.decoherence_C_param, prms.decoherence_eps_param, prms.isNBRA);       
   }
 
   else if(prms.decoherence_times_type==2){
@@ -720,14 +743,13 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
   if(prms.dephasing_informed==1){
     Eadi = get_Eadi(ham); 
     MATRIX ave_gaps(*prms.ave_gaps);
-    dephasing_informed_correction(decoherence_rates, Eadi, ave_gaps);
+    dephasing_informed_correction(decoherence_rates, Eadi, ave_gaps, prms.isNBRA);
   }
 
-
   //============ Apply decoherence corrections ==================
   // SDM and alike methods
   if(prms.decoherence_algo==0){
-    Coeff = sdm(Coeff, prms.dt, act_states, decoherence_rates, prms.sdm_norm_tolerance);
+    Coeff = sdm(Coeff, prms.dt, act_states, decoherence_rates, prms.sdm_norm_tolerance, prms.isNBRA);
   }
   // BCSH
   else if(prms.decoherence_algo==3){ 
@@ -737,7 +759,7 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
 
   else if(prms.decoherence_algo==4){
 //    MATRIX decoh_rates(ndof, ntraj);
-    Coeff = mfsd(p, Coeff, invM, prms.dt, decoherence_rates, ham, rnd);
+    Coeff = mfsd(p, Coeff, invM, prms.dt, decoherence_rates, ham, rnd, prms.isNBRA);
   }
 
 //  exit(0);
@@ -801,6 +823,7 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
 
     /// New version, as of 8/3/2020
     vector<int> old_states(act_states);
+    //cout << "Flag before dish" << endl;
     act_states = dish(prms, q, p, invM, Coeff, projectors, ham, act_states, coherence_time, decoherence_rates, rnd);
 
     /// Velocity rescaling
@@ -841,7 +864,6 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
   }
 
 //    exit(0);
-
   St.clear();
   Eadi.clear();
   decoherence_rates.clear();
@@ -852,13 +874,13 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMA
 
 
 //    exit(0);
-
 }
 
 
 
 
 
+
 }// namespace libdyn
 }// liblibra
 

src/dyn/Dynamics.h

@@ -74,8 +74,11 @@ void update_Hamiltonian_p(bp::dict prms, nHamiltonian& ham, MATRIX& p, MATRIX& i
 
 CMATRIX transform_amplitudes(int rep_in, int rep_out, CMATRIX& C, nHamiltonian& ham);
 
-
+// Adding the NBRA flag to the functions in the header
+vector<CMATRIX> compute_St(nHamiltonian& ham, int isNBRA);
 vector<CMATRIX> compute_St(nHamiltonian& ham);
+
+vector<CMATRIX> compute_St(nHamiltonian& ham, vector<CMATRIX>& Uprev, int isNBRA);
 vector<CMATRIX> compute_St(nHamiltonian& ham, vector<CMATRIX>& Uprev);
 
 

src/dyn/dyn_control_params.cpp

@@ -53,7 +53,7 @@ dyn_control_params::dyn_control_params(){
   convergence = 0;
   max_number_attempts = 100;
   min_probability_reordering = 0.0;
-
+  isNBRA = 0;
 
   ///================= Surface hopping: proposal, acceptance =======================
   tsh_method = 0;
@@ -187,6 +187,7 @@ void dyn_control_params::set_parameters(bp::dict params){
     else if(key=="convergence") { convergence = bp::extract<int>(params.values()[i]);  }
     else if(key=="max_number_attempts") { max_number_attempts = bp::extract<int>(params.values()[i]);  }
     else if(key=="min_probability_reordering") { min_probability_reordering = bp::extract<double>(params.values()[i]);  }
+    else if(key=="isNBRA") { isNBRA = bp::extract<int>(params.values()[i]);   }
 
     ///================= Surface hopping: proposal, acceptance =======================
     else if(key=="tsh_method") { tsh_method = bp::extract<int>(params.values()[i]);  }