diff --git a/CMakeLists.txt b/CMakeLists.txt index 95862058d..efdaede88 100755 --- a/CMakeLists.txt +++ b/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}) diff --git a/README.md b/README.md index 5fc3ae3a2..25ad39714 100755 --- a/README.md +++ b/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 diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index c0d4e5a83..62f7d9538 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -248,7 +248,7 @@ CMATRIX transform_amplitudes(int rep_in, int rep_out, CMATRIX& C, nHamiltonian& //vector compute_St(nHamiltonian& ham, CMATRIX** Uprev){ -vector compute_St(nHamiltonian& ham, vector& Uprev){ +vector compute_St(nHamiltonian& ham, vector& Uprev, int isNBRA){ /** This function computes the time-overlap matrices for all trajectories @@ -258,16 +258,26 @@ vector compute_St(nHamiltonian& ham, vector& Uprev){ int ntraj = ham.children.size(); vector St(ntraj, CMATRIX(nst, nst)); - + if(isNBRA==1){ + St[0] = Uprev[0].H() * ham.children[0]->get_basis_transform(); + } + else{ for(int traj=0; trajget_basis_transform(); } - + } return St; } -vector compute_St(nHamiltonian& ham){ +vector compute_St(nHamiltonian& ham, vector& Uprev){ + + int is_nbra = 0; + return compute_St(ham, Uprev, is_nbra); + +} + +vector compute_St(nHamiltonian& ham, int isNBRA){ /** This function computes the time-overlap matrices for all trajectories @@ -277,16 +287,23 @@ vector compute_St(nHamiltonian& ham){ int ntraj = ham.children.size(); vector St(ntraj, CMATRIX(nst, nst)); - + if(isNBRA==1){ + St[0] = ham.children[0]->get_time_overlap_adi(); + } + else{ for(int traj=0; trajget_time_overlap_adi(); } - + } return St; } +vector compute_St(nHamiltonian& ham){ + int is_nbra = 1; + return compute_St(ham, is_nbra); +} void apply_afssh(dyn_variables& dyn_var, CMATRIX& C, vector& act_states, MATRIX& invM, @@ -487,7 +504,6 @@ void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector project_out_states(ntraj); // for DISH vector Uprev; - vector St(ntraj, CMATRIX(nst, nst)); - vector Eadi(ntraj, CMATRIX(nst, nst)); - vector decoherence_rates(ntraj, MATRIX(nst, nst)); - vector Ekin(ntraj, 0.0); - vector 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 St(ntraj1, CMATRIX(nst, nst)); + vector Eadi(ntraj1, CMATRIX(nst, nst)); + vector decoherence_rates(ntraj1, MATRIX(nst, nst)); + vector Ekin(ntraj1, 0.0); + vector prev_ham_dia(ntraj1, MATRIX(nst, nst)); MATRIX gamma(ndof, ntraj); MATRIX p_traj(ndof, 1); vector t1(ndof, 0); for(dof=0;dofget_ham_dia().real(); + //} + //else{ + for(traj=0; trajget_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 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, vectornadi, 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 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 compute_St(nHamiltonian& ham, int isNBRA); vector compute_St(nHamiltonian& ham); + +vector compute_St(nHamiltonian& ham, vector& Uprev, int isNBRA); vector compute_St(nHamiltonian& ham, vector& Uprev); diff --git a/src/dyn/dyn_control_params.cpp b/src/dyn/dyn_control_params.cpp index 589671e16..22465610b 100755 --- a/src/dyn/dyn_control_params.cpp +++ b/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(params.values()[i]); } else if(key=="max_number_attempts") { max_number_attempts = bp::extract(params.values()[i]); } else if(key=="min_probability_reordering") { min_probability_reordering = bp::extract(params.values()[i]); } + else if(key=="isNBRA") { isNBRA = bp::extract(params.values()[i]); } ///================= Surface hopping: proposal, acceptance ======================= else if(key=="tsh_method") { tsh_method = bp::extract(params.values()[i]); } diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index 432b666f2..70a48f040 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -220,6 +220,7 @@ class dyn_control_params{ */ double min_probability_reordering; + ///=============================================================================== ///================= Surface hopping: proposal, acceptance ======================= ///=============================================================================== @@ -412,6 +413,18 @@ class dyn_control_params{ MATRIX* ave_gaps; + /** + A flag for NBRA calculations. Since in NBRA, the Hamiltonian is the same for all the trajectories + we can only compute the Hamiltonian related properties once for one trajectory and increase the speed of calculations. + If we set the value to 1 it will consider the NBRA type calculations and other integers the Hamiltonian related properties + are computed for all trajectories. + + Options: + - 0: no NBRA - Hamiltonians for all trajectories are computed explicitly [ default ] + - 1: the NBRA is involved - the calculations of the Hamiltonian are conducted for only 1 trajectory, + and re-used by all other SH trajectories. + */ + int isNBRA; ///=============================================================================== diff --git a/src/dyn/dyn_decoherence.h b/src/dyn/dyn_decoherence.h index 4ebf85658..b523924ce 100755 --- a/src/dyn/dyn_decoherence.h +++ b/src/dyn/dyn_decoherence.h @@ -39,9 +39,11 @@ namespace libdyn{ CMATRIX sdm(CMATRIX& Coeff, double dt, int act_st, MATRIX& decoh_rates, double tol); CMATRIX sdm(CMATRIX& Coeff, double dt, int act_st, MATRIX& decoh_rates); +CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates, double tol, int isNBRA); CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates, double tol); CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates); + void project_out(CMATRIX& Coeff, int traj, int i); void collapse(CMATRIX& Coeff, int traj, int i, int collapse_option); @@ -61,17 +63,27 @@ CMATRIX bcsh(CMATRIX& Coeff, double dt, vector& act_states, MATRIX& reversa // For MF-SD of Schwartz +CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& decoherence_rates, nHamiltonian& ham, Random& rnd, int isNBRA); CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& decoherence_rates, nHamiltonian& ham, Random& rnd); ///================ In dyn_decoherence_time.cpp =================================== +MATRIX edc_rates(CMATRIX& Hvib, double Ekin, double C_param, double eps_param, int isNBRA); MATRIX edc_rates(CMATRIX& Hvib, double Ekin, double C_param, double eps_param); + +vector edc_rates(vector& Hvib, vector& Ekin, double C_param, double eps_param, int isNBRA); vector edc_rates(vector& Hvib, vector& Ekin, double C_param, double eps_param); + +void dephasing_informed_correction(MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps, int isNBRA); +void dephasing_informed_correction(vector& decoh_rates, vector& Hvib, MATRIX& ave_gaps, int isNBRA); + void dephasing_informed_correction(MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps); void dephasing_informed_correction(vector& decoh_rates, vector& Hvib, MATRIX& ave_gaps); + + MATRIX coherence_intervals(CMATRIX& Coeff, MATRIX& rates); MATRIX coherence_intervals(CMATRIX& Coeff, vector& rates); diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 9dc8dd3dc..203ad4b49 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -37,8 +37,6 @@ CMATRIX sdm(CMATRIX& Coeff, double dt, int act_st, MATRIX& decoh_rates, double t \param[in] act_st [ integer ] The active state index \param[in] decoh_rates [ MATRIX ] The matrix of decoherence (pure dephasing) rates between all pairs of states \param[in] tol [double] The maximal acceptable deviation of the p_aa_old from 1. If the p_aa_old < 1.0 + tol, then renormalize it to 1.0 - - The function returns: C [ CMATRIX ] - the updated state of the electronic DOF, in the same data type as the input */ @@ -152,7 +150,8 @@ CMATRIX sdm(CMATRIX& Coeff, double dt, int act_st, MATRIX& decoh_rates){ } -CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates, double tol){ + +CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates, double tol, int isNBRA){ /** \brief The generic framework of the Simplified Decay of Mixing (SDM) method of Granucci, G.; Persico, M. J. Chem. Phys. 2007, 126, 134114) @@ -164,6 +163,7 @@ CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& deco \param[in] act_st [ integer ] The active state index \param[in] decoh_rates [ MATRIX ] The matrix of decoherence (pure dephasing) rates between all pairs of states \param[in] tol [double] The maximal acceptable deviation of the p_aa_old from 1. If the p_aa_old < 1.0 + tol, then renormalize it to 1.0 + \param[in] isNBRA [integer] If this flag is set to 1, then the Hamiltonian related properties are only computed for one of the trajectories. The function returns: # C [ CMATRIX ] - the updated state of the electronic DOF, in the same data type as the input @@ -180,26 +180,37 @@ CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& deco CMATRIX coeff(nadi, 1); CMATRIX res(nadi, ntraj); - - for(traj=0; traj& act_st, vector& decoh_rates, double tol){ + int is_nbra = 0; + return sdm(Coeff, dt, act_st, decoh_rates, tol, is_nbra); + +} + CMATRIX sdm(CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates){ - - double tol = 0.0; - return sdm(Coeff, dt, act_st, decoh_rates, tol); + double tol = 0.0; + int is_nbra = 0; + return sdm(Coeff, dt, act_st, decoh_rates, tol, is_nbra); + } + void project_out(CMATRIX& Coeff, int traj, int i){ /** Projects the state `i` out of a coherent superposition of states @@ -640,7 +651,7 @@ CMATRIX bcsh(CMATRIX& Coeff, double dt, vector& act_states, MATRIX& reversa -CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& decoherence_rates, nHamiltonian& ham, Random& rnd){ +CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& decoherence_rates, nHamiltonian& ham, Random& rnd, int isNBRA){ /** \brief Mean field with stochastic decoherence @@ -682,8 +693,14 @@ CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& // Probability of decoherence on state i // Here, we assume that the state-only decoherence rates are on the diagonal - double P_i = (std::conj(c_i) * c_i).real() * decoherence_rates[itraj].get(i, i) * dt; - + double P_i; + if(isNBRA==1){ + P_i = (std::conj(c_i) * c_i).real() * decoherence_rates[0].get(i, i) * dt; + } + else + { + P_i = (std::conj(c_i) * c_i).real() * decoherence_rates[itraj].get(i, i) * dt; + } if(ksi& +CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& decoherence_rates, nHamiltonian& ham, Random& rnd){ + + int is_nbra = 0; + return mfsd(p, Coeff, invM, dt, decoherence_rates, ham, rnd, is_nbra); + +} + + }// namespace libdyn diff --git a/src/dyn/dyn_decoherence_time.cpp b/src/dyn/dyn_decoherence_time.cpp index 328a29f52..44887d33e 100755 --- a/src/dyn/dyn_decoherence_time.cpp +++ b/src/dyn/dyn_decoherence_time.cpp @@ -24,7 +24,7 @@ namespace liblibra{ namespace libdyn{ -MATRIX edc_rates(CMATRIX& Hvib, double Ekin, double C_param, double eps_param){ +MATRIX edc_rates(CMATRIX& Hvib, double Ekin, double C_param, double eps_param, int isNBRA){ /** This function computes the decoherence rates matrix used in the energy-based decoherence scheme of Granucci-Persico and Truhlar @@ -34,7 +34,7 @@ MATRIX edc_rates(CMATRIX& Hvib, double Ekin, double C_param, double eps_param){ \param[in] Ekin [ float ] The classical kinetic energy of nuclei. Units = Ha \param[in] C_param [ float ] The method parameter, typically set to 1.0 Ha \param[in] eps_param [ float ] The method parameter, typically set to 0.1 Ha - + \param[in] isNBRA [ int ] The method for considering NBRA-type calculations */ int i,j; @@ -55,11 +55,29 @@ MATRIX edc_rates(CMATRIX& Hvib, double Ekin, double C_param, double eps_param){ } -vector edc_rates(vector& Hvib, vector& Ekin, double C_param, double eps_param){ +MATRIX edc_rates(CMATRIX& Hvib, double Ekin, double C_param, double eps_param){ + int is_nbra = 0; + return edc_rates(Hvib, Ekin, C_param, eps_param, is_nbra); +} + + + +vector edc_rates(vector& Hvib, vector& Ekin, double C_param, double eps_param, int isNBRA){ int ntraj = Hvib.size(); int nst = Hvib[0].n_cols; + // A new variable instead of ntraj + int ntraj1; + + if(isNBRA==1){ + ntraj1 = 1; + } + else{ + ntraj1 = Hvib.size(); + } + vector res(ntraj1, MATRIX(nst, nst)); + if(isNBRA==1){ if(Ekin.size()!=ntraj){ cout<<"ERROR in edc_rates: the sizes of the input variables Hvib and Ekin are inconsistent\n"; cout<<"Hvib.size() = "< edc_rates(vector& Hvib, vector& Ekin, double C_p cout<<"exiting...\n"; exit(0); } - - vector res(ntraj, MATRIX(nst, nst)); - for(int traj=0; traj edc_rates(vector& Hvib, vector& Ekin, double C_param, double eps_param){ + int is_nbra = 0; + return edc_rates(Hvib, Ekin, C_param, eps_param, is_nbra); +} -void dephasing_informed_correction(MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps){ + + + +void dephasing_informed_correction(MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps, int isNBRA){ /** This function computes the corrected dephasing rates The correction is based on the instnataneous energy levels and on the average @@ -89,6 +112,7 @@ void dephasing_informed_correction(MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& a \param[in, out] decoh_rates [ MATRIX ] uncorrected decoherence rates [units: a.u.t.^-1] \param[in] Hvib [ CMATRIX ] Instantaneous vibronic Hamiltonian [units: Ha] \param[in] ave_gaps [ MATRIX ] time-averaged module of the energy level gaps: <|E_i - E_j|> [units: Ha] + \param[in] isNBRA [ int ] The method for considering NBRA-type calculations */ @@ -117,10 +141,22 @@ void dephasing_informed_correction(MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& a } -void dephasing_informed_correction(vector& decoh_rates, vector& Hvib, MATRIX& ave_gaps){ +void dephasing_informed_correction(MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps){ + + int is_nbra = 0; + dephasing_informed_correction(decoh_rates, Hvib, ave_gaps, is_nbra); +} + + + +void dephasing_informed_correction(vector& decoh_rates, vector& Hvib, MATRIX& ave_gaps, int isNBRA){ int ntraj = Hvib.size(); + if(isNBRA==1){ + dephasing_informed_correction(decoh_rates[0], Hvib[0], ave_gaps, isNBRA); + } + else{ if(decoh_rates.size()!=ntraj){ cout<<"ERROR in dephasing_informed_correction: the sizes of the input variables \ decoh_rates and Hvib are inconsistent\n"; @@ -132,9 +168,17 @@ void dephasing_informed_correction(vector& decoh_rates, vector& for(int traj=0; traj& decoh_rates, vector& Hvib, MATRIX& ave_gaps){ + + int is_nbra = 0; + dephasing_informed_correction(decoh_rates, Hvib, ave_gaps, is_nbra); } diff --git a/src/dyn/dyn_hop_acceptance.cpp b/src/dyn/dyn_hop_acceptance.cpp index 42fa7dece..35fb59bad 100755 --- a/src/dyn/dyn_hop_acceptance.cpp +++ b/src/dyn/dyn_hop_acceptance.cpp @@ -445,6 +445,7 @@ vector accept_hops(dyn_control_params& prms, int ntraj_active = which_trajectories.size(); int nst = C.n_rows; int itraj, traj, dof, i, idof; + int isNBRA = prms.isNBRA; vector fstates(ntraj, -1); MATRIX p_tr(ndof, 1); @@ -476,8 +477,18 @@ vector accept_hops(dyn_control_params& prms, p_tr = p.col(traj); double T_i = compute_kinetic_energy(p_tr, invM); // initial kinetic energy - hvib = ham.children[traj]->get_ham_adi(); - hvib = projectors[traj].H() * hvib * projectors[traj]; + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(itraj==0){ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } + } + + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy @@ -540,10 +551,17 @@ vector accept_hops(dyn_control_params& prms, int new_st = proposed_states[traj]; if(old_st != new_st){ - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(itraj==0){ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; - + } + } + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy @@ -596,9 +614,17 @@ vector accept_hops(dyn_control_params& prms, int new_st = proposed_states[traj]; if(old_st != new_st){ - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(itraj==0){ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; + } + } + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy @@ -636,10 +662,17 @@ vector accept_hops(dyn_control_params& prms, int old_st = initial_states[traj]; int new_st = proposed_states[traj]; - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(itraj==0){ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; - + } + } + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy @@ -664,10 +697,17 @@ vector accept_hops(dyn_control_params& prms, int old_st = initial_states[traj]; int new_st = proposed_states[traj]; - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(itraj==0){ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; - + } + } + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy @@ -693,10 +733,17 @@ vector accept_hops(dyn_control_params& prms, int old_st = initial_states[traj]; int new_st = proposed_states[traj]; - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(itraj==0){ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; - + } + } + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy @@ -816,7 +863,7 @@ void handle_hops_nuclear(dyn_control_params& prms, int ntraj = q.n_cols; int nst = C.n_rows; int traj, idof, dof, i; - + int isNBRA = prms.isNBRA; MATRIX p_tr(ndof, 1); CMATRIX hvib(nst, nst); CMATRIX nac(nst, nst); @@ -837,10 +884,17 @@ void handle_hops_nuclear(dyn_control_params& prms, p_tr = p.col(traj); double T_i = compute_kinetic_energy(p_tr, invM, which_dofs); // initial kinetic energy - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(traj==0){ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; - + } + } + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy double T_f = T_i + E_i - E_f; // predicted final kinetic energy @@ -907,10 +961,17 @@ void handle_hops_nuclear(dyn_control_params& prms, int old_st = old_states[traj]; int new_st = new_states[traj]; - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(traj==0){ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; - + } + } + else{ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy @@ -960,9 +1021,17 @@ void handle_hops_nuclear(dyn_control_params& prms, int old_st = old_states[traj]; int new_st = new_states[traj]; - + // Do the calculations only for itraj=0 if isNBRA is 1 + if(isNBRA==1){ + if(traj==0){ + hvib = ham.children[traj]->get_ham_adi(); + hvib = projectors[traj].H() * hvib * projectors[traj]; + } + } + else{ hvib = ham.children[traj]->get_ham_adi(); hvib = projectors[traj].H() * hvib * projectors[traj]; + } double E_i = hvib.get(old_st, old_st).real(); // initial potential energy double E_f = hvib.get(new_st, new_st).real(); // final potential energy diff --git a/src/dyn/dyn_hop_proposal.cpp b/src/dyn/dyn_hop_proposal.cpp index 17a65a4d8..8cb8be129 100755 --- a/src/dyn/dyn_hop_proposal.cpp +++ b/src/dyn/dyn_hop_proposal.cpp @@ -327,6 +327,7 @@ vector hop_proposal_probabilities(dyn_control_params& prms, int ntraj = q.n_cols; int nst = C.n_rows; int traj, dof, i; + int isNBRA = prms.isNBRA; vector nucl_stenc_x(ndof, 0); for(i=0;i hop_proposal_probabilities(dyn_control_params& prms, pop_submatrix(C, coeff, el_stenc_x, el_stenc_y); + if(isNBRA==1){ + // Only compute the Hvib for one traj, traj==0 + if(traj==0){ + Hvib = ham.children[traj]->get_hvib_adi(); + + //Transform Hamiltonian to the dynamically-consistent form: + Hvib = projectors[traj].H() * Hvib * projectors[traj]; + } + } + + else{ + // Compute the Hvib for all traj Hvib = ham.children[traj]->get_hvib_adi(); //Transform Hamiltonian to the dynamically-consistent form: Hvib = projectors[traj].H() * Hvib * projectors[traj]; - + } if(prms.tsh_method == 0){ // FSSH diff --git a/src/dyn/dyn_projectors.cpp b/src/dyn/dyn_projectors.cpp index 17e14bccc..9a4a2ed5e 100755 --- a/src/dyn/dyn_projectors.cpp +++ b/src/dyn/dyn_projectors.cpp @@ -568,6 +568,10 @@ void update_projectors(dyn_control_params& prms, vector& projectors, CMATRIX st(nst, nst); CMATRIX ist(nst, nst); CMATRIX projector_old(nst, nst); + // Instead of setting many of-else in the for loop we can set ntraj=1 + if(prms.isNBRA==1){ + ntraj = 1; + } for(int traj=0; traj& ham, int rep); void propagate_electronic(double dt, CMATRIX& C, vector& projector, vector& ham, int rep); +void propagate_electronic(double dt, CMATRIX& C, vector& projector, vector& ham, int rep, int isNBRA); + void propagate_electronic(double dt, CMATRIX& C, nHamiltonian& ham, int rep, int level); void propagate_electronic(double dt, CMATRIX& C, vector& projector, nHamiltonian& ham, int rep, int level); //void propagate_electronic(double dt, nHamiltonian& ham, int rep); diff --git a/src/dyn/electronic/Electronic_Dynamics1.cpp b/src/dyn/electronic/Electronic_Dynamics1.cpp index ba8d5ac61..afbb285da 100755 --- a/src/dyn/electronic/Electronic_Dynamics1.cpp +++ b/src/dyn/electronic/Electronic_Dynamics1.cpp @@ -137,7 +137,6 @@ void Electronic::propagate_electronic(double dt,Hamiltonian* ham){ \param[in] ham The pointer to Hamiltonian object, that affects the dynamics */ - libdyn::libelectronic::propagate_electronic(dt,this,ham); } @@ -156,7 +155,6 @@ void Electronic::propagate_electronic(double dt,Hamiltonian& ham){ \param[in] ham The reference to Hamiltonian object, that affects the dynamics */ - libdyn::libelectronic::propagate_electronic(dt,this,&ham); } @@ -178,7 +176,6 @@ void propagate_electronic(double dt,Electronic* el,Hamiltonian* ham){ */ - int i,j; double dt_half = 0.5*dt; @@ -908,7 +905,7 @@ void propagate_electronic(double dt, CMATRIX& Coeff, CMATRIX& Hvib, CMATRIX& S){ void propagate_electronic(double dt, CMATRIX& C, nHamiltonian& ham, int rep){ - if(rep==0){ // diabatic + if(rep==0){ // diabatic CMATRIX Hvib(ham.ndia, ham.ndia); Hvib = ham.get_hvib_dia(); CMATRIX Sdia(ham.ndia, ham.ndia); Sdia = ham.get_ovlp_dia(); @@ -1032,15 +1029,45 @@ void propagate_electronic(double dt, CMATRIX& C, vector& ham, int } -void propagate_electronic(double dt, CMATRIX& C, vector& projector, vector& ham, int rep){ +void propagate_electronic(double dt, CMATRIX& C, vector& projector, vector& ham, int rep, int isNBRA){ + if(isNBRA!=1){ if(C.n_cols!=ham.size()){ cout<<"ERROR in void propagate_electronic(double dt, CMATRIX& C, vector& ham, int rep): \n"; cout<<"C.n_cols = "<nadi, ham[0]->nadi); Hvib = ham[0]->get_hvib_adi(); + int sz = Hvib.n_cols; + CMATRIX* I; I = new CMATRIX(sz, sz); I->load_identity(); + CMATRIX* C1; C1 = new CMATRIX(sz, sz); *C1 = complex(0.0, 0.0); // eigenvectors + CMATRIX* Heig; Heig = new CMATRIX(sz, sz); *Heig = complex(0.0,0.0); // eigenvalues + CMATRIX* expH; expH = new CMATRIX(sz, sz); *expH = complex(0.0,0.0); + + libmeigen::solve_eigen(Hvib, *I, *Heig, *C1, 0); // Hvib_eff * C1 = I * C1 * Heig ==> ham[0] = C1 * Heig * C.H() + complex one(0.0, 1.0); + for(i=0;i val = std::exp(-one*Heig->get(i,i)*dt ); + expH->set(i,i,val); + } + + *expH = (*C1) * (*expH) * ((*C1).H()); + // Instead of the for loop + C = *expH * C; + + delete expH; + delete I; + delete C1; + delete Heig; + + } + else { int nst = C.n_rows; int ntraj = C.n_cols; @@ -1054,10 +1081,16 @@ void propagate_electronic(double dt, CMATRIX& C, vector& projector, vec for(int st=0; st& projector, vector& ham, int rep){ + + int is_nbra = 0; + propagate_electronic(dt, C, projector, ham, rep); + +} void propagate_electronic(double dt, CMATRIX& C, nHamiltonian& ham, int rep, int level){ diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index 0e2713018..50b7485a2 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -71,6 +71,7 @@ void export_dyn_control_params_objects(){ .def_readwrite("MK_verbosity", &dyn_control_params::MK_verbosity) .def_readwrite("convergence", &dyn_control_params::convergence) .def_readwrite("max_number_attempts", &dyn_control_params::max_number_attempts) + .def_readwrite("isNBRA", &dyn_control_params::isNBRA) ///================= Surface hopping: proposal, acceptance ======================= .def_readwrite("tsh_method", &dyn_control_params::tsh_method) @@ -161,14 +162,20 @@ void export_dyn_decoherence_objects(){ (CMATRIX& Coeff, double dt, int act_st, MATRIX& decoh_rates) = &sdm; def("sdm", expt_sdm_v2); + CMATRIX (*expt_sdm_v3) - (CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates, double tol) = &sdm; + (CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates, double tol, int isNBRA) = &sdm; def("sdm", expt_sdm_v3); CMATRIX (*expt_sdm_v4) - (CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates) = &sdm; + (CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates, double tol) = &sdm; def("sdm", expt_sdm_v4); + CMATRIX (*expt_sdm_v5) + (CMATRIX& Coeff, double dt, vector& act_st, vector& decoh_rates) = &sdm; + def("sdm", expt_sdm_v5); + + void (*expt_project_out_v1)(CMATRIX& Coeff, int traj, int i) = &project_out; @@ -195,31 +202,58 @@ void export_dyn_decoherence_objects(){ CMATRIX (*expt_mfsd_v1) + (MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& decoherence_rates, + nHamiltonian& ham, Random& rnd, int isNBRA) = &mfsd; + def("mfsd", expt_mfsd_v1); + + CMATRIX (*expt_mfsd_v2) (MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& decoherence_rates, nHamiltonian& ham, Random& rnd) = &mfsd; + def("mfsd", expt_mfsd_v2); ///================ In dyn_decoherence_time.cpp =================================== MATRIX (*expt_edc_rates_v1) - (CMATRIX& Hvib, double Ekin, double C_param, double eps_param) = &edc_rates; + (CMATRIX& Hvib, double Ekin, double C_param, double eps_param, int isNBRA) = &edc_rates; def("edc_rates", expt_edc_rates_v1); - vector (*expt_edc_rates_v2) - (vector& Hvib, vector& Ekin, - double C_param, double eps_param) = &edc_rates; + MATRIX (*expt_edc_rates_v2) + (CMATRIX& Hvib, double Ekin, double C_param, double eps_param) = &edc_rates; def("edc_rates", expt_edc_rates_v2); + vector (*expt_edc_rates_v3) + (vector& Hvib, vector& Ekin, + double C_param, double eps_param, int isNBRA) = &edc_rates; + def("edc_rates", expt_edc_rates_v3); + + vector (*expt_edc_rates_v4) + (vector& Hvib, vector& Ekin, + double C_param, double eps_param, int isNBRA) = &edc_rates; + def("edc_rates", expt_edc_rates_v4); + + + void (*expt_dephasing_informed_correction_v1) - (MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps) = &dephasing_informed_correction; + (MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps, int isNBRA) = &dephasing_informed_correction; def("dephasing_informed_correction", expt_dephasing_informed_correction_v1); void (*expt_dephasing_informed_correction_v2) - (vector& decoh_rates, vector& Hvib, MATRIX& ave_gaps) = &dephasing_informed_correction; + (MATRIX& decoh_rates, CMATRIX& Hvib, MATRIX& ave_gaps) = &dephasing_informed_correction; def("dephasing_informed_correction", expt_dephasing_informed_correction_v2); + + void (*expt_dephasing_informed_correction_v3) + (vector& decoh_rates, vector& Hvib, MATRIX& ave_gaps, int isNBRA) = &dephasing_informed_correction; + def("dephasing_informed_correction", expt_dephasing_informed_correction_v3); + + void (*expt_dephasing_informed_correction_v4) + (vector& decoh_rates, vector& Hvib, MATRIX& ave_gaps) = &dephasing_informed_correction; + def("dephasing_informed_correction", expt_dephasing_informed_correction_v4); + + MATRIX (*expt_coherence_intervals_v1)(CMATRIX& Coeff, MATRIX& rates) = &coherence_intervals; def("coherence_intervals", expt_coherence_intervals_v1); @@ -645,10 +679,15 @@ void export_Dyn_objects(){ - vector (*expt_compute_St_v1)(nHamiltonian& ham) = &compute_St; + vector (*expt_compute_St_v1)(nHamiltonian& ham, int isNBRA) = &compute_St; def("compute_St", expt_compute_St_v1); - vector (*expt_compute_St_v2)(nHamiltonian& ham, vector& Uprev) = &compute_St; + vector (*expt_compute_St_v2)(nHamiltonian& ham) = &compute_St; def("compute_St", expt_compute_St_v2); + vector (*expt_compute_St_v3)(nHamiltonian& ham, vector& Uprev, int isNBRA) = &compute_St; + def("compute_St", expt_compute_St_v3); + vector (*expt_compute_St_v4)(nHamiltonian& ham, vector& Uprev) = &compute_St; + def("compute_St", expt_compute_St_v4); + diff --git a/src/dyn/qtag/qtag.cpp b/src/dyn/qtag/qtag.cpp index 847f84483..96f6450f8 100755 --- a/src/dyn/qtag/qtag.cpp +++ b/src/dyn/qtag/qtag.cpp @@ -349,6 +349,188 @@ complex LHA(CMATRIX* Ham1, CMATRIX* Ham2, }// LHA +complex LHAe(int i, int j, + MATRIX& q1, MATRIX& p1, MATRIX& s1, MATRIX& alp1, int n1, + MATRIX& q2, MATRIX& p2, MATRIX& s2, MATRIX& alp2, int n2, + nHamiltonian& ham){ +/** + """Returns the (complex) value for the potential *v* on an energetic surface specified by *nsurf* from two basis + functions defined by their parameters *qpasi* and *qpasj*, respectively. The computation employs the Local Harmonic + Approximation (LHA), which requires a potential function *pot* as well as its first and second derivatives. + The potential parameters are stored in the dict 'model_params' as defined in qtag_config. + + Args: + univ (dictionary): Dictionary containing various system parameters. + + nstate (integer): Integer specifying the potential surface to be calculated (0 = ground, 1 = first excited, ...) + + qpasi (list): The ndof-by-4 parameter list of the i-th basis function. Each entry is an ndof-by-1 column MATRIX. + + qpasj (list): The ndof-by-4 parameter list of the j-th basis function. Each entry is an ndof-by-1 column MATRIX. + + params (dictionary): Dictionary containing the potential parameters. + + libra_model (function object): Function object containing the Libra potential model for the ground and excited states. + + Returns: + v (complex): The complex potential v computed via the local harmonic approximation between basis functions i and j. + """ +*/ + + complex v; + int ndof = q1.n_rows; + +//Holstein coupling + double A = 1.0; + double B = 1.5811; + double C = 2.0; + +//Tully1 coupling +// double A = 0.005; +// double B = 1.0; +// double C = 0.0; + +//Tully2 coupling +// double A = 0.015; +// double B = 0.06; +// double C = 0.0; + + if(n1==n2){// LHA for single-surface elements + + v = LHA(ham.children[i]->ham_dia, ham.children[j]->ham_dia, + ham.children[i]->d1ham_dia, ham.children[j]->d1ham_dia, + ham.children[i]->d2ham_dia, ham.children[j]->d2ham_dia, + q1, p1, s1, alp1, n1, q2, p2, s2, alp2, n2); + + } + else if(n1!=n2){// Exact integral for Gaussian coupling + + v = (0.0, 0.0); + + for(int dof=0; dof expt(aCq1*aCq1+aCq2*aCq2-dp*dp+2.0*aCq1*aCq2, 2.0*dp*(aCq1+aCq2)); + v += prefac1*exp(prefac2*expt); + }// for dof + }//end if + + return v; + +}//LHAe + +complex BATe(int i, int j, + MATRIX& q1, MATRIX& p1, MATRIX& s1, MATRIX& alp1, int n1, + MATRIX& q2, MATRIX& p2, MATRIX& s2, MATRIX& alp2, int n2, + nHamiltonian& ham){ +/** + """Returns the (complex) value for the potential *v* on an energetic surface specified by *nsurf* from two basis + functions defined by their parameters *qpasi* and *qpasj*, respectively. The computation employs the Local Harmonic + Approximation (LHA), which requires a potential function *pot* as well as its first and second derivatives. + The potential parameters are stored in the dict 'model_params' as defined in qtag_config. + + Args: + univ (dictionary): Dictionary containing various system parameters. + + nstate (integer): Integer specifying the potential surface to be calculated (0 = ground, 1 = first excited, ...) + + qpasi (list): The ndof-by-4 parameter list of the i-th basis function. Each entry is an ndof-by-1 column MATRIX. + + qpasj (list): The ndof-by-4 parameter list of the j-th basis function. Each entry is an ndof-by-1 column MATRIX. + + params (dictionary): Dictionary containing the potential parameters. + + libra_model (function object): Function object containing the Libra potential model for the ground and excited states. + + Returns: + v (complex): The complex potential v computed via the local harmonic approximation between basis functions i and j. + """ +*/ + + complex v; + int ndof = q1.n_rows; + +//Holstein coupling + double A = 1.0; + double B = 1.5811; + double C = 2.0; + +//Tully1 coupling +// double A = 0.005; +// double B = 1.0; +// double C = 0.0; + +//Tully2 coupling +// double A = 0.015; +// double B = 0.06; +// double C = 0.0; + + if(n1==n2){// BAT for single-surface elements + + v = BAT(ham.children[i]->ham_dia, ham.children[j]->ham_dia, + ham.children[i]->d1ham_dia, ham.children[j]->d1ham_dia, + q1, p1, s1, alp1, n1, q2, p2, s2, alp2, n2); + + } + else if(n1!=n2){// Exact integral for Gaussian coupling + + v = (0.0, 0.0); + + for(int dof=0; dof expt(aCq1*aCq1+aCq2*aCq2-dp*dp+2.0*aCq1*aCq2, 2.0*dp*(aCq1+aCq2)); + v += prefac1*exp(prefac2*expt); + }// for dof + }//end if + + return v; + +}//BATe CMATRIX qtag_potential(MATRIX& q1, MATRIX& p1, MATRIX& s1, MATRIX& alp1, int n1, vector& traj_on_surf_n1, MATRIX& q2, MATRIX& p2, MATRIX& s2, MATRIX& alp2, int n2, vector& traj_on_surf_n2, @@ -402,8 +584,14 @@ CMATRIX qtag_potential(MATRIX& q1, MATRIX& p1, MATRIX& s1, MATRIX& alp1, int n1, ham.children[i]->d2ham_dia, ham.children[j]->d2ham_dia, qi, pi, si, ai, n1, qj, pj, sj, aj, n2); }// LHA - - //cout<<"i= "<=nfiles: + index -= nfiles + + model_params.update({"timestep":index}) + #model_params.update({"timestep":i}) if rep_tdse==0: compute_dynamics(q, p, iM, Cdia, projectors, states, ham, compute_model, model_params, dyn_params, rnd, therm, dyn_var) elif rep_tdse==1: compute_dynamics(q, p, iM, Cadi, projectors, states, ham, compute_model, model_params, dyn_params, rnd, therm, dyn_var) -# sys.exit(0) - if _savers["txt_saver"]!=None: _savers["txt_saver"].save_data_txt( F"{prefix}", properties_to_save, "a", i) @@ -1087,489 +1126,6 @@ def run_dynamics(_q, _p, _iM, _Cdia, _Cadi, _projectors, _states, _dyn_params, c - -def run_dynamics_old(_q, _p, _iM, _Cdia, _Cadi, _projectors, _states, _dyn_params, compute_model, _model_params, rnd): - """ - - This is a deprecated function - it returns/saves the results of the calculations via either - text files or memory, both ways are pretty bad. But if for some strange reason you can't use HDF5 - you may want to use this version. So, I keep int here, just in case. - - This function is similar to the one above, but it stores the computed properties in files - - Args: - _q ( MATRIX(nnucl, ntraj) ): coordinates of the "classical" particles [units: Bohr] - _p ( MATRIX(nnucl, ntraj) ): momenta of the "classical" particles [units: a.u. of momenta] - _iM ( MATRIX(nnucl, 1) ): masses of classical particles [units: a.u.^-1] - _Cdia ( CMATRIX(ndia, ntraj) ): amplitudes of the diabatic basis states - _Cadi ( CMATRIX(nadi, ntraj) ): amplitudes of the adiabatic basis states - _projectors ( list of CMATRIX(nadi, nadi) ): cumulative phase correction and state tracking matrices - _states ( intList, or list of ntraj ints ): the quantum state of each trajectory - - dyn_params ( dictionary ): parameters controlling the execution of the dynamics - Can contain: - - * **dyn_params["rep_tdse"]** ( int ): selects the representation in which - nuclear/electronic (Ehrenfest core) dynamics is executed. The representation - used to integrate TD-SE - - - 0: diabatic representation - - 1: adiabatic representation [ default ] - - - * **dyn_params["rep_ham"]** (int): The representation of the Hamiltonian update: - - - 0: diabatic [ default ] - - 1: adiabatic - - This is the representation in which the computed properties are assumed to be - For instance, we may have set it to 1, to read the adiabatic energies and couplings, - to bypass the diabatic-to-adiabatic transformation, which may be useful in some atomistic - calculations, or with the NBRA - - - * **dyn_params["rep_sh"]** ( int ): selects the representation which is - used to perform surface hopping - - - 0: diabatic representation - - 1: adiabatic representation [ default ] - - - * **dyn_params["rep_lz"]** ( int ): The representation to compute LZ probabilitieis: - - - 0: diabatic [ default ] - - 1: adiabatic - - - * **dyn_params["tsh_method"]** ( int ): Formula for computing SH probabilities: - - - -1: adiabatic - no hops [ default ] - - 0: FSSH - - 1: GFSH - - 2: MSSH - - 3: DISH - - - * **dyn_params["force_method"]** ( int ): How to compute forces in the dynamics: - - - 0: don't compute forces at all - e.g. we do not really need them - - 1: state-specific as in the TSH or adiabatic (including adiabatic excited states) [ default ] - - 2: Ehrenfest - - - * **dyn_params["nac_update_method"]** ( int ): How to update NACs and vibronic Hamiltonian before - electronic TD-SE propagation: - - - 0: don't update them (e.g. for simplest NAC) - - 1: update according to changed momentum and existing derivative couplings [ default ] - - - * **dyn_params["rep_force"]** ( int ): In which representation to compute forces: - - - 0: diabatic - - 1: adiabatic [ default ] - - - * **dyn_params["use_boltz_factor"]** ( int ): Whether to scale the SH probabilities - by the Boltzmann factor: - - - 0: do not scale [ default ] - - 1: scale - - - * **dyn_params["Temperature"]** ( double ): Temperature of the system [ default: 300 K ] - - - * **dyn_params["do_reverse"]** ( int ): what to do with velocities on frustrated hops: - - - 0: do not revert momenta at the frustrated hops - - 1: do revert the momenta [ default ] - - * **dyn_params["vel_rescale_opt"]** ( int ): How to rescale momenta if the hops are successful: - - - -1: do not rescale, as in the NBRA [ default ] - - 0: rescale in the diabatic basis - don't care about the - velocity directions, just a uniform rescaling - - 1: rescale along the directions of derivative couplings - - - * **dyn_params["do_phase_correction"]** ( int ): the algorithm to correct phases on adiabatic states - - - 0: no phase correction - - 1: according to our phase correction algorithm [ default: 1 ] - - - * **dyn_params["tol"]** ( double ): the minimal value of the time-overlap for considering - phase corrections - no correction applied if the time-overlap is smaller in magnitude than - this parameter [ default: 1e-3 ] - - - * **dyn_params["state_tracking_algo"]** ( int ): the algorithm to keep track of the states' identities - - - 0: no state tracking - - 1: Sato - - 2: using the mincost, Munkres-Kuhn [ default: 2 ] - - - * **dyn_params["MK_alpha"]** ( double ): Munkres-Kuhn alpha - (selects the range of orbitals included in reordering) [default: 0.0] - - - * **dyn_params["MK_verbosity"]** ( double ): Munkres-Kuhn verbosity: - - - 0: no extra output [ default ] - - 1: print details - - - - * **dyn_params["entanglement_opt"]** ( int ): A selector of a method to couple the trajectories in this ensemble: - - - 0: no coupling across trajectories [ default ] - - 1: ETHD - - 2: ETHD3 (experimental) - - 22: another flavor of ETHD3 (experimental) - - * **dyn_params["ETHD3_alpha"]** ( double ): Gaussian exponents that dresses up the trajectories in the ETHD3 method - in the coordinate space, that is ~exp(-alpha*(R-R0)^2 ) [ default: 0.0 ] - - - * **dyn_params["ETHD3_beta"]** ( double ): Gaussian exponents that dresses up the trajectories in the ETHD3 method - in the coordinate space, that is ~exp(-alpha*(P-P0)^2 ) [ default: 0.0 ] - - - * **dyn_params["decoherence_algo"]** ( int ): selector of the method to incorporate decoherence: - - - -1: no decoherence [ default ] - - 0: MSDM - - 1: ID-A - - - * **dyn_params["decoh_rates"]** ( MATRIX(ntates, nstates) ): the matrix of dephasing rates [ units: a.u. of time ^-1 ] - - - * **dyn_params["dt"]** ( double ): the nuclear and electronic integration - timestep [ units: a.u. of time, default: 41.0 ] - - - * **dyn_params["nsteps"]** ( int ): the number of NA-MD steps to do [ default: 1 ] - - - * **dyn_params["prefix"]** ( string ): the name of the folder to be created by this function. - All the data files will be created in that folder - - - * **dyn_params["output_level"]** ( int ): controls what info to return as the result of this function - - - -1: all return values are None [ default ] - - 0: also, all return values are none - - 1: 1-D info - energies, SE and SH populations averaged over ensemble vs. time - - 2: 2-D info - coordinates, momenta, SE amplitudes, and SH states populations - for each individual trajectory vs. time - - 3: 3-D info - St, Hvib_adi, Hvib_dia matrices (energies, couplings, etc.) for each - individual trajectory vs. time - - - * **dyn_params["file_output_level"]** ( int ): controls what info to print out into files - - - -1: print nothing at all [ default ] - - 0: 0-D info - just the parameters of the simulation - - 1: 1-D info - energies, SE and SH populations averaged over ensemble vs. time - - 2: 2-D info - coordinates, momenta, SE amplitudes, and SH states populations - for each individual trajectory vs. time - - 3: 3-D info - St, Hvib_adi, Hvib_dia matrices (energies, couplings, etc.) for each - individual trajectory vs. time - - - compute_model ( PyObject ): the pointer to the Python function that performs the Hamiltonian calculations - - _model_params ( dictionary ): contains the selection of a model and the parameters - for that model Hamiltonian - - rnd ( Random ): random numbers generator object - - - - Returns: - tuple: with the elements of the tuple listed below. - - * obs_T ( list of `nsteps` doubles ): time [units: a.u.] - * obs_q ( list of `nsteps` MATRIX(nnucl, ntraj) ): coordinates of all trajectories [ units: Bohr ] - * obs_p ( list of `nsteps` MATRIX(nnucl, ntraj) ): momenta of all trajectories [ units: a.u. ] - * obs_Ekin ( list of `nsteps` doubles ): average kinetic energy of an ensemble of trajectories [units: a.u.] - * obs_Epot ( list of `nsteps` doubles ): average potential energy of an ensemble of trajectories [units: a.u.] - * obs_Etot ( list of `nsteps` doubles ): average total energy of an ensemble of trajectories [units: a.u.] - * obs_dEkin ( list of `nsteps` doubles ): standard deviation of kinetic energy of an ensemble of trajectories [units: a.u.] - * obs_dEpot ( list of `nsteps` doubles ): standard deviation of potential energy of an ensemble of trajectories [units: a.u.] - * obs_dEtot ( list of `nsteps` doubles ): standard deviation of total energy of an ensemble of trajectories [units: a.u.] - * obs_Cadi ( list of `nsteps` CMATRIX(nadi, ntraj) ): amplitudes of adiabatic electronic states of all trajectories - * obs_Cdia ( list of `nsteps` CMATRIX(ndia, ntraj) ): amplitudes of diabatic electronic states of all trajectories - * obs_dm_adi ( list of `nsteps` CMATRIX(nadi, nadi) ): ensemble-averaged density matrix in adiabatic basis - * obs_dm_dia ( list of `nsteps` CMATRIX(ndia, ndia) ): ensemble-averaged density matrix in diabatic basis - * obs_pop ( list of `nsteps` MATRIX(nadi, 1) ): ensemble-averaged TSH populations of adiabatic states - * obs_states ( list of `nsteps` of lists of `ntraj` ints): # indices of the quantum states of each trajectory - * obs_hvib_adi ( list of `ntraj` lists, each being a list of `nsteps` objects CMATRIX(nadi, nadi) ): trajectory-resolved - vibronic Hamiltonians for each timestep in the adiabatic representation - * obs_hvib_dia ( list of `ntraj` lists, each being a list of `nsteps` objects CMATRIX(ndia, ndia) ): trajectory-resolved - vibronic Hamiltonians for each timestep in the diabatic representation - * obs_St ( list of `ntraj` lists, each being a list of `nsteps` objects CMATRIX(nadi, nadi) ): trajectory-resolved - time-overlaps of the adiabatic wavefunctions - * obs_U ( list of `ntraj` lists, each being a list of `nsteps` objects CMATRIX(ndia, nadi) ): trajectory-resolved - diabatic-to-adiabatic transformation matrix - - - .. note:: - the elements are None, if they are excluded by the input optinos - - """ - - - # Create copies of the input dynamical variables, so we could run several such - # functions with the same input variables without worries that they will be altered - # inside of each other - - model_params = dict(_model_params) - dyn_params = dict(_dyn_params) - - q = MATRIX(_q) - p = MATRIX(_p) - iM = MATRIX(_iM) - Cdia = CMATRIX(_Cdia) - Cadi = CMATRIX(_Cadi) - states = intList() - projectors = CMATRIXList() - - - for i in range(len(_states)): - states.append(_states[i]) - projectors.append(CMATRIX(_projectors[i])) - - DR = MATRIX(len(_states), len(_states)) - AG = MATRIX(len(_states), len(_states)) - - - # Parameters and dimensions - critical_params = [ ] - default_params = { "rep_tdse":1, "rep_ham":0, "rep_sh":1, "rep_lz":0, "tsh_method":-1, - "force_method":1, "nac_update_method":1, "rep_force":1, - "use_boltz_factor":0, "Temperature":300.0, "do_reverse":1, "vel_rescale_opt":-1, - "do_phase_correction":1, "tol":1e-3, - "state_tracking_algo":2, "MK_alpha":0.0, "MK_verbosity":0, - "entanglement_opt":0, "ETHD3_alpha":0.0, "ETHD3_beta":0.0, - "decoherence_algo":-1, "decoherence_rates":DR, - "decoherence_times_type":0, "decoherence_C_param":1.0, - "decoherence_eps_param":0.1, "dephasing_informed":0, - "ave_gaps":AG, "instantaneous_decoherence_variant":1, "collapse_option":0, - "ensemble":0, "thermostat_params":{}, - "dt":1.0*units.fs2au, "nsteps":1, - "output_level":-1, "file_output_level":-1, "prefix":"tmp" - } - - comn.check_input(dyn_params, default_params, critical_params) - - prefix = dyn_params["prefix"] - rep_tdse = dyn_params["rep_tdse"] - nsteps = dyn_params["nsteps"] - dt = dyn_params["dt"] - tol = dyn_params["tol"] - output_level = dyn_params["output_level"] - file_output_level = dyn_params["file_output_level"] - do_phase_correction = dyn_params["do_phase_correction"] - state_tracking_algo = dyn_params["state_tracking_algo"] - force_method = dyn_params["force_method"] - - ndia = Cdia.num_of_rows - nadi = Cadi.num_of_rows - nnucl= q.num_of_rows - ntraj= q.num_of_cols - - - - obs_T = None - obs_q, obs_p = None, None - obs_Ekin, obs_Epot, obs_Etot = None, None, None - obs_dEkin, obs_dEpot, obs_dEtot = None, None, None - obs_Cadi, obs_Cdia = None, None - obs_dm_adi, obs_dm_dia = None, None - obs_pop, obs_states = None, None - obs_hvib_adi, obs_hvib_dia, obs_St, obs_U = None, None, None, None - - if output_level>=1: - obs_T = [] # time - obs_Ekin = [] # average kinetic energy - obs_Epot = [] # average potential energy - obs_Etot = [] # average total energy - obs_dEkin = [] # kinetic energy fluctuation - obs_dEpot = [] # potential energy fluctuation - obs_dEtot = [] # total energy fluctuation - obs_dm_adi = [] # average SE-based density matrix in adiabatic basis - obs_dm_dia = [] # average SE-based density matrix in diabatic basis - obs_pop = [] # average SH-based populations adiabatic basis - - if output_level>=2: - obs_q = [] # coordinates of all trajectories - obs_p = [] # momenta of all trajectories - obs_Cadi = [] # average TD-SE amplitudes in the adiabatic basis - obs_Cdia = [] # average TD-SE amplitudes in the diabatic basis - obs_states = [] # indices of the quantum states of each trajectory - - if output_level>=3: - obs_hvib_adi = [] # vibronic Hamiltonians for each trajectory in adiabatic rep. - obs_hvib_dia = [] # vibronic Hamiltonians for each trajectory in diabatic rep. - obs_St = [] # time-overlaps for each trajectory in the adiabatic rep. - obs_U = [] # dia-to-adia transformation matrices for each trajectory - - for tr in range(ntraj): - obs_hvib_adi.append([]) - obs_hvib_dia.append([]) - obs_St.append([]) - obs_U.append([]) - - - # Create an output directory, if not present - if file_output_level>=0: - if not os.path.isdir(prefix): - os.mkdir(prefix) - - # Simulation parameters - f = open("%s/_dyn_params.txt" % (prefix),"w") - f.write( str(dyn_params) ); f.close() - - f = open("%s/_model_params.txt" % (prefix),"w") - f.write( str(model_params) ); f.close() - - # Ensemble-resolved info (so 1D) - if file_output_level >= 1: - f = open("%s/energies.txt" % (prefix), "w"); f.close() # average kinetic, potential, and total energies and their fluctuations - f = open("%s/D_adi.txt" % (prefix), "w"); f.close() # average SE-based density matrix in adiabatic basis - f = open("%s/D_dia.txt" % (prefix), "w"); f.close() # average SE-based density matrix in diabatic basis - f = open("%s/SH_pop.txt" % (prefix), "w"); f.close() # average SH-based populations adiabatic basis - - # Trajectory-resolved 1D info (so 2D) - if file_output_level >= 2: - f = open("%s/q.txt" % (prefix), "w"); f.close() # coordinates of all trajectories - f = open("%s/p.txt" % (prefix), "w"); f.close() # momenta of all trajectories - f = open("%s/C_adi.txt" % (prefix), "w"); f.close() # TD-SE amplitudes in the adiabatic basis - f = open("%s/C_dia.txt" % (prefix), "w"); f.close() # TD-SE amplitudes in the diabatic basis - f = open("%s/states.txt" % (prefix), "w"); f.close() # indices of the quantum states of each trajectory - - if file_output_level >= 3: - # Trajectory-resolved 2D info (so 3D) - for tr in range(ntraj): - f = open("%s/Hvib_adi_%i.txt" % (prefix, tr), "w"); f.close() # vibronic Hamiltonians for each trajectory in adiabatic rep. - f = open("%s/Hvib_dia_%i.txt" % (prefix, tr), "w"); f.close() # vibronic Hamiltonians for each trajectory in diabatic rep. - f = open("%s/St_%i.txt" % (prefix, tr), "w"); f.close() # time-overlaps along each trajectory - f = open("%s/basis_transform_%i.txt" % (prefix, tr), "w"); f.close() # dia-to-adi transformationfor each trajectory. - - - # ======= Hierarchy of Hamiltonians ======= - ham = nHamiltonian(ndia, nadi, nnucl) - ham.add_new_children(ndia, nadi, nnucl, ntraj) - ham.init_all(2,1) - model_params.update({"timestep":0}) - - update_Hamiltonian_q(dyn_params, q, projectors, ham, compute_model, model_params) - update_Hamiltonian_p(dyn_params, ham, p, iM) - - - U = [] - for tr in range(ntraj): - U.append(ham.get_basis_transform(Py2Cpp_int([0, tr]) )) - - - - # Do the propagation - for i in range(nsteps): - - #============ Compute and output properties =========== - # Amplitudes, Density matrix, and Populations - if rep_tdse==0: - # Diabatic to raw adiabatic - ham.ampl_dia2adi(Cdia, Cadi, 0, 1) - # Raw adiabatic to dynamically-consistent adiabatic - Cadi = dynconsyst_to_raw(Cadi, projectors) - - elif rep_tdse==1: - ham.ampl_adi2dia(Cdia, Cadi, 0, 1) - - dm_dia, dm_adi = tsh_stat.compute_dm(ham, Cdia, Cadi, projectors, rep_tdse, 1) - pops = tsh_stat.compute_sh_statistics(nadi, states) - - - # Energies - if force_method in [0, 1]: - Ekin, Epot, Etot, dEkin, dEpot, dEtot = tsh_stat.compute_etot_tsh(ham, p, Cdia, Cadi, projectors, states, iM, rep_tdse) - elif force_method in [2]: - Ekin, Epot, Etot, dEkin, dEpot, dEtot = tsh_stat.compute_etot(ham, p, Cdia, Cadi, projectors, iM, rep_tdse) - - - # ============== Debug ==================== - nrm = (Cadi.H() * Cadi ).tr() - - x = 0.0 - for traj in range(ntraj): - x += (projectors[traj].H() * projectors[traj]).tr() - - #print(F"step = {i} nrm = {nrm} projectors measure = {x}") - #============ End of Debug ================= - - # Memory output - if output_level >= 1: - obs_T.append(i*dt) - obs_Ekin.append(Ekin) - obs_Epot.append(Epot) - obs_Etot.append(Etot) - obs_dEkin.append(dEkin) - obs_dEpot.append(dEpot) - obs_dEtot.append(dEtot) - obs_dm_adi.append(CMATRIX(dm_adi)) - obs_dm_dia.append(CMATRIX(dm_dia)) - obs_pop.append(MATRIX(pops)) - - # Memory output - if output_level >= 2: - obs_q.append(MATRIX(q)) - obs_p.append(MATRIX(p)) - obs_Cadi.append(CMATRIX(Cadi)) - obs_Cdia.append(CMATRIX(Cdia)) - obs_states.append(list(states)) - - # File output - res12 = q, p, Ekin, Epot, Etot, dEkin, dEpot, dEtot, Cadi, Cdia, dm_adi, dm_dia, pops, states - dynamics_io.print_results12(i, dt, res12, prefix, file_output_level) - - - for tr in range(ntraj): - x = ham.get_basis_transform(Py2Cpp_int([0, tr]) ) - St = U[tr] * x.H() - U[tr] = CMATRIX(x) - - # Memory output - if output_level >= 3: - obs_hvib_adi[tr].append( CMATRIX(ham.get_hvib_adi(Py2Cpp_int([0, tr])) ) ) - obs_hvib_dia[tr].append( CMATRIX(ham.get_hvib_dia(Py2Cpp_int([0, tr])) ) ) - obs_St[tr].append( CMATRIX(St) ) - obs_U[tr].append( CMATRIX(x) ) - - # File output - res3 = ham.get_hvib_adi(Py2Cpp_int([0, tr])), ham.get_hvib_dia(Py2Cpp_int([0, tr])), St, U[tr] - dynamics_io.print_results3(i, dt, res3, prefix, file_output_level, tr) - - - #============ Propagate =========== - model_params.update({"timestep":i}) - if rep_tdse==0: - compute_dynamics(q, p, iM, Cdia, projectors, states, ham, compute_model, model_params, dyn_params, rnd) - elif rep_tdse==1: - compute_dynamics(q, p, iM, Cadi, projectors, states, ham, compute_model, model_params, dyn_params, rnd) - - - - return obs_T, obs_q, obs_p, obs_Ekin, obs_Epot, obs_Etot, obs_dEkin, obs_dEpot, obs_dEtot, \ - obs_Cadi, obs_Cdia, obs_dm_adi, obs_dm_dia, obs_pop, obs_states, obs_hvib_adi, obs_hvib_dia, obs_St, obs_U - - - - def generic_recipe(q, p, iM, _dyn_params, compute_model, _model_params, _init_elec, rnd): """ This function initializes electronic variables and Hamiltonians for a given set of @@ -1607,16 +1163,23 @@ def generic_recipe(q, p, iM, _dyn_params, compute_model, _model_params, _init_el """ - comn.check_input(_model_params, { }, [ "model0" ] ) - comn.check_input(_dyn_params, { "rep_tdse":1 }, [ ] ) + model_params = dict(_model_params) + dyn_params = dict(_dyn_params) + init_elec = dict(_init_elec) + + comn.check_input( model_params, { }, [ "model0" ] ) + comn.check_input( dyn_params, { "rep_tdse":1, "is_nbra":0 }, [ ] ) + is_nbra = dyn_params["is_nbra"] + + init_elec.update({"is_nbra":is_nbra}) # Internal parameters nnucl, ntraj = q.num_of_rows, q.num_of_cols # Initialize electronic variables - either diabatic or adiabatic - Cdia, Cadi, projectors, states = init_electronic_dyn_var(_init_elec, rnd) + Cdia, Cadi, projectors, states = init_electronic_dyn_var(init_elec, rnd) ndia, nadi = Cdia.num_of_rows, Cadi.num_of_rows @@ -1624,26 +1187,32 @@ def generic_recipe(q, p, iM, _dyn_params, compute_model, _model_params, _init_el # compute the diabatic-to-adiabatic transformation matrices and # transform the amplitudes accordingly ham = nHamiltonian(ndia, nadi, nnucl) - ham.add_new_children(ndia, nadi, nnucl, ntraj) + if is_nbra == 1: + ham.add_new_children(ndia, nadi, nnucl, 1) + else: + ham.add_new_children(ndia, nadi, nnucl, ntraj) ham.init_all(2,1) - if _init_elec["rep"]==0: - if _dyn_params["rep_tdse"]==1: - model_params = dict(_model_params) - model_params.update({"model":_model_params["model0"]}) - update_Hamiltonian_q({"rep_tdse":1, "rep_ham":0}, q, projectors, ham, compute_model, model_params ) + if init_elec["rep"]==0: + if dyn_params["rep_tdse"]==1: + model_params1 = dict(model_params) + model_params1.update({"model":model_params["model0"]}) + update_Hamiltonian_q({"rep_tdse":1, "rep_ham":0}, q, projectors, ham, compute_model, model_params1 ) Cadi = transform_amplitudes(0, 1, Cdia, ham) - elif _init_elec["rep"]==1: - if _dyn_params["rep_tdse"]==0: - model_params = dict(_model_params) - model_params.update({"model":_model_params["model0"]}) - update_Hamiltonian_q({"rep_tdse":1, "rep_ham":0}, q, projectors, ham, compute_model, model_params ) + elif init_elec["rep"]==1: + if dyn_params["rep_tdse"]==0: + model_params1 = dict(model_params) + model_params1.update({"model":model_params["model0"]}) + update_Hamiltonian_q({"rep_tdse":1, "rep_ham":0}, q, projectors, ham, compute_model, model_params1 ) Cdia = transform_amplitudes(1, 0, Cdia, ham) - res = run_dynamics(q, p, iM, Cdia, Cadi, projectors, states, _dyn_params, compute_model, _model_params, rnd) + #if _dyn_params["isNBRA"]==1: + # res = run_dynamics_nbra(q, p, iM, Cdia, Cadi, projectors, states, _dyn_params, compute_model, _model_params, rnd) + #else: + res = run_dynamics(q, p, iM, Cdia, Cadi, projectors, states, _dyn_params, compute_model, model_params, rnd) return res diff --git a/src/libra_py/models/Holstein.py b/src/libra_py/models/Holstein.py index fbf0b1a26..3f6b1e973 100755 --- a/src/libra_py/models/Holstein.py +++ b/src/libra_py/models/Holstein.py @@ -413,6 +413,8 @@ def Holstein5(q, params, full_id): * obj.ovlp_dia ( CMATRIX(4,4) ): overlap of the basis (diabatic) states [ identity ] * obj.d1ham_dia ( list of 1 CMATRIX(4,4) objects ): derivatives of the diabatic Hamiltonian w.r.t. the nuclear coordinate + * obj.d2ham_dia ( list of 1 CMATRIX(4,4) objects ): + second derivatives of the diabatic Hamiltonian w.r.t. the nuclear coordinate * obj.dc1_dia ( list of 1 CMATRIX(4,4) objects ): derivative coupling in the diabatic basis [ zero ] """ @@ -439,6 +441,7 @@ def Holstein5(q, params, full_id): Hdia = CMATRIX(n,n) Sdia = CMATRIX(n,n) d1ham_dia = CMATRIXList(); d1ham_dia.append( CMATRIX(n,n) ) + d2ham_dia = CMATRIXList(); d2ham_dia.append( CMATRIX(n,n) ) dc1_dia = CMATRIXList(); dc1_dia.append( CMATRIX(n,n) ) Id = Cpp2Py(full_id) @@ -467,12 +470,23 @@ def Holstein5(q, params, full_id): if i!=j: d1ham_dia[k].set(i,j, -2.0*alpha[i][j] * (x-x_nm[i][j]) * V[i][j] * math.exp(-alpha[i][j] * (x-x_nm[i][j])**2 ) * (1.0+0.0j) ) + for k in [0]: + # d2 Hdia / dR_0 2 + for i in range(n): + d2ham_dia[k].set(i,i, k_n[i]*(1.0+0.0j) ) + + for k in [0]: + for i in range(n): + for j in range(n): + if i!=j: + d2ham_dia[k].set(i,j, -2.0*alpha[i][j] * V[i][j] * math.exp(-alpha[i][j] * (x-x_nm[i][j])**2 )*(1.0-2.0*alpha[i][j]*(x-x_nm[i][j])**2)*(1.0+0.0j) ) obj = tmp() obj.ham_dia = Hdia obj.ovlp_dia = Sdia obj.d1ham_dia = d1ham_dia + obj.d2ham_dia = d2ham_dia obj.dc1_dia = dc1_dia return obj diff --git a/src/libra_py/models/Zhu.py b/src/libra_py/models/Zhu.py index 2a4f20308..6a5c9b2f7 100755 --- a/src/libra_py/models/Zhu.py +++ b/src/libra_py/models/Zhu.py @@ -32,7 +32,7 @@ class tmp: pass -def dual_RZD(q, params): +def dual_RZD(q, params, full_id): """ Dual Rosen-Zener-Demkov model @@ -81,8 +81,12 @@ def dual_RZD(q, params): Sdia = CMATRIX(2,2) d1ham_dia = CMATRIXList(); d1ham_dia.append( CMATRIX(2,2) ) dc1_dia = CMATRIXList(); dc1_dia.append( CMATRIX(2,2) ) - - x = q.get(0) + + Id = Cpp2Py(full_id) + indx = Id[-1] + x = q.get(0, indx) + + Sdia.set(0,0, 1.0+0.0j); Sdia.set(0,1, 0.0+0.0j); Sdia.set(1,0, 0.0+0.0j); Sdia.set(1,1, 1.0+0.0j); @@ -119,7 +123,7 @@ def dual_RZD(q, params): -def dual_LZS(q, params): +def dual_LZS(q, params, full_id): """ Dual Landau-Zener-Stuckelberg model @@ -170,8 +174,11 @@ def dual_LZS(q, params): Sdia = CMATRIX(2,2) d1ham_dia = CMATRIXList(); d1ham_dia.append( CMATRIX(2,2) ) dc1_dia = CMATRIXList(); dc1_dia.append( CMATRIX(2,2) ) + + Id = Cpp2Py(full_id) + indx = Id[-1] + x = q.get(0, indx) - x = q.get(0) Sdia.set(0,0, 1.0+0.0j); Sdia.set(0,1, 0.0+0.0j); Sdia.set(1,0, 0.0+0.0j); Sdia.set(1,1, 1.0+0.0j); @@ -211,7 +218,7 @@ def dual_LZS(q, params): -def Renner_Teller(q, params): +def Renner_Teller(q, params, full_id): """ Renner-Teller @@ -257,8 +264,11 @@ def Renner_Teller(q, params): Sdia = CMATRIX(2,2) d1ham_dia = CMATRIXList(); d1ham_dia.append( CMATRIX(2,2) ) dc1_dia = CMATRIXList(); dc1_dia.append( CMATRIX(2,2) ) - - x = q.get(0) + + Id = Cpp2Py(full_id) + indx = Id[-1] + x = q.get(0, indx) + Sdia.set(0,0, 1.0+0.0j); Sdia.set(0,1, 0.0+0.0j); Sdia.set(1,0, 0.0+0.0j); Sdia.set(1,1, 1.0+0.0j); diff --git a/src/libra_py/tsh_stat.py b/src/libra_py/tsh_stat.py index 2f2974116..8eef70fb4 100755 --- a/src/libra_py/tsh_stat.py +++ b/src/libra_py/tsh_stat.py @@ -32,6 +32,7 @@ import sys import math import copy +import numpy as np if sys.platform=="cygwin": from cyglibra_core import * @@ -42,7 +43,7 @@ from . import units from . import probabilities - +from libra_py import data_conv def compute_etot(ham, p, Cdia, Cadi, projectors, iM, rep): """Computes the Ehrenfest potential energy @@ -221,7 +222,7 @@ def compute_etot_tsh(ham, p, Cdia, Cadi, projectors, act_states, iM, rep): -def compute_dm(ham, Cdia, Cadi, projectors, rep, lvl): +def compute_dm(ham, Cdia, Cadi, projectors, rep, lvl, isNBRA=0): """ Compute the trajectory-averaged density matrices in diabatic @@ -245,6 +246,10 @@ def compute_dm(ham, Cdia, Cadi, projectors, rep, lvl): - 0: ham is the actual Hamiltonian to use (use with single trajectory), - 1: ham is the parent of the Hamiltonians to use (use with multiple trajectories) + isNBRA ( int ): The flag for NBRA type calculations: + - 0: the Hamiltonian related properties are computed for all of the trajectories [default] + - 1: the Hamiltonian related properties are computed only for one trajectory + Returns: tuple: ( dm_dia, dm_adi ): @@ -267,8 +272,10 @@ def compute_dm(ham, Cdia, Cadi, projectors, rep, lvl): # Raw dm_dia_raw, dm_adi_raw = CMATRIX(ndia, ndia), CMATRIX(nadi, nadi) + if isNBRA==1: + + traj = 0 - for traj in range(0,ntraj): indx = None if lvl==0: indx = Py2Cpp_int([0]) @@ -276,10 +283,10 @@ def compute_dm(ham, Cdia, Cadi, projectors, rep, lvl): indx = Py2Cpp_int([0,traj]) S = ham.get_ovlp_dia(indx) - U = ham.get_basis_transform(indx) - + U = ham.get_basis_transform(indx) + if rep==0: - + dm_tmp = S * Cdia.col(traj) * Cdia.col(traj).H() * S # Dia dyn-consistent @@ -293,27 +300,75 @@ def compute_dm(ham, Cdia, Cadi, projectors, rep, lvl): # Adi raw dm_adi_raw = dm_adi_raw + projectors[traj] * tmp * projectors[traj].H() - - + + elif rep==1: su = S * U # Raw - c = Cadi.col(traj) - dm_tmp = c * c.H() - + #c = Cadi.col(traj) + dm_tmp = Cadi.col(traj) * Cadi.col(traj).H() + # Adi dynamically-consistent DM dm_adi = dm_adi + dm_tmp - + # Adi raw dm_tmp_raw = projectors[traj] * dm_tmp * projectors[traj].H() dm_adi_raw = dm_adi_raw + dm_tmp_raw # Dia dynamically-consistent DM - dm_dia = dm_dia + su * dm_tmp * su.H() + dm_dia = dm_dia + su * dm_tmp * su.H() # Dia raw - dm_dia_raw = dm_dia_raw + su * dm_tmp_raw * su.H() + dm_dia_raw = dm_dia_raw + su * dm_tmp_raw * su.H() + + else: + for traj in range(0,ntraj): + indx = None + if lvl==0: + indx = Py2Cpp_int([0]) + elif lvl==1: + indx = Py2Cpp_int([0,traj]) + + S = ham.get_ovlp_dia(indx) + U = ham.get_basis_transform(indx) + + if rep==0: + + dm_tmp = S * Cdia.col(traj) * Cdia.col(traj).H() * S + + # Dia dyn-consistent + dm_dia = dm_dia + dm_tmp + + # Dia raw - i'm not sure about that one, so lets keep it just zero for now + + # Adi dyn-consistent + tmp = U.H() * dm_tmp * U + dm_adi = dm_adi + tmp + + # Adi raw + dm_adi_raw = dm_adi_raw + projectors[traj] * tmp * projectors[traj].H() + + + elif rep==1: + su = S * U + + # Raw + c = Cadi.col(traj) + dm_tmp = c * c.H() + + # Adi dynamically-consistent DM + dm_adi = dm_adi + dm_tmp + + # Adi raw + dm_tmp_raw = projectors[traj] * dm_tmp * projectors[traj].H() + dm_adi_raw = dm_adi_raw + dm_tmp_raw + + # Dia dynamically-consistent DM + dm_dia = dm_dia + su * dm_tmp * su.H() + + # Dia raw + dm_dia_raw = dm_dia_raw + su * dm_tmp_raw * su.H() dm_dia = dm_dia / float(ntraj) @@ -324,8 +379,7 @@ def compute_dm(ham, Cdia, Cadi, projectors, rep, lvl): return dm_dia, dm_adi, dm_dia_raw, dm_adi_raw - -def compute_sh_statistics(nstates, istate, projectors): +def compute_sh_statistics(nstates, istate, projectors, isNBRA=0): """ This function computes the SH statistics for an ensemble of trajectories @@ -338,6 +392,9 @@ def compute_sh_statistics(nstates, istate, projectors): Each element of the list is the state index for that trajectory. In other words, istate[0] is the quantum state for a trajectory 0, istate[1] is the quantum state for a trajectory 1, etc. + isNBRA ( int ): The flag for NBRA type calculations: + - 0: the Hamiltonian related properties are computed for all of the trajectories [ default ] + - 1: the Hamiltonian related properties are computed only for one trajectory Returns: MATRIX(nstates, 1): coeff_sh: The list containing the average @@ -348,20 +405,34 @@ def compute_sh_statistics(nstates, istate, projectors): ntraj = len(istate) coeff_sh = MATRIX(nstates, 1) coeff_sh_raw = MATRIX(nstates, 1) - - for i in range(0,ntraj): - st = istate[i] - pop = CMATRIX(nstates, nstates) - pop.set(st, st, 1.0+0.0j) - - pop_raw = projectors[i] * pop * projectors[i].H() - - for j in range(nstates): - coeff_sh.add(j, 0, pop.get(j,j).real) - coeff_sh_raw.add(j, 0, pop_raw.get(j,j).real) + + if isNBRA==1: + + pops_numpy = np.zeros((nstates, ntraj)) + for i in range(ntraj): + pops_numpy[istate[i], i] = 1 + + coeff_sh_numpy = np.average(pops_numpy, axis=1).reshape(nstates,1) + coeff_sh = data_conv.nparray2MATRIX(coeff_sh_numpy) + coeff_sh_raw = coeff_sh + + else: + + for i in range(0,ntraj): + + st = istate[i] + #print(st) + pop = CMATRIX(nstates, nstates) + pop.set(st, st, 1.0+0.0j) + #projectors[i].show_matrix() + pop_raw = projectors[i] * pop * projectors[i].H() - coeff_sh.scale(-1,0, 1.0/float(ntraj)) - coeff_sh_raw.scale(-1,0, 1.0/float(ntraj)) + for j in range(nstates): + coeff_sh.add(j, 0, pop.get(j,j).real) + coeff_sh_raw.add(j, 0, pop_raw.get(j,j).real) + + coeff_sh.scale(-1,0, 1.0/float(ntraj)) + coeff_sh_raw.scale(-1,0, 1.0/float(ntraj)) return coeff_sh, coeff_sh_raw diff --git a/src/libra_py/workflows/nbra/step2.py b/src/libra_py/workflows/nbra/step2.py index efe0c16cc..1dfecbd9e 100755 --- a/src/libra_py/workflows/nbra/step2.py +++ b/src/libra_py/workflows/nbra/step2.py @@ -408,33 +408,44 @@ def run_cp2k_libint_step2(params): os.system(F"mkdir {params['all_logfiles']}") if not os.path.exists(params['all_pdosfiles']): os.system(F"mkdir {params['all_pdosfiles']}") + # setting up the initial step and final step istep = params['istep'] fstep = params['fstep'] + # spherical or cartesian GTOs is_spherical = params['is_spherical'] + # number of processors nprocs = params['nprocs'] + # CP2K executable cp2k_exe = params['cp2k_exe'] + # mpirun executable mpi_executable = params['mpi_executable'] + # Unrestricted spin calculations isUKS = params['isUKS'] + # Extended tight-binding calculations isxTB = params['isxTB'] + # Periodic calculation falg is_periodic = params['is_periodic'] + # Cube visualization cube_visualization = params['cube_visualization'] + if cube_visualization: - # Read the tcl file data + # Read the tcl file data states_to_plot = params['states_to_plot'] if isUKS: phase_factors_alpha = np.ones(( len(states_to_plot), 1 )) phase_factors_beta = np.ones(( len(states_to_plot), 1 )) else: phase_factors_alpha = np.ones(( len(states_to_plot), 1 )) + # the counter for a job steps, this counter is needed for not reading # the data of a molden file twice counter = 0 @@ -455,6 +466,7 @@ def run_cp2k_libint_step2(params): os.system(F'{mpi_executable} -n {nprocs} {cp2k_exe} -i Diag_libra-{step}.inp -o step_{step}.log') molden_filename = F'Diag_libra-{step}-1_0.molden' print('Done with step', step,'Elapsed time:',time.time()-t1) + # now if the counter is equal to zero # just compute the MO overlap of that step. if counter == 0: @@ -483,6 +495,7 @@ def run_cp2k_libint_step2(params): print(F'Computing the AO overlaps between R({translational_vector[0]},{translational_vector[1]},{translational_vector[2]}) and R(0,0,0)') shell_1p, l_vals = molden_methods.molden_file_to_libint_shell(molden_filename, is_spherical, is_periodic, cell, translational_vector) AO_S += compute_overlaps(shell_1,shell_1p, nprocs) + print('Done with computing atomic orbital overlaps. Elapsed time:', time.time()-t1) t1 = time.time() print('Turning the MATRIX to numpy array...') @@ -496,6 +509,7 @@ def run_cp2k_libint_step2(params): t1 = time.time() new_indices = CP2K_methods.resort_molog_eigenvectors(l_vals) eigenvectors_1 = [] + for j in range(len(eig_vect_1)): # the new and sorted eigenvector eigenvector_1 = eig_vect_1[j] @@ -580,6 +594,7 @@ def run_cp2k_libint_step2(params): shell_2p, l_vals = molden_methods.molden_file_to_libint_shell(molden_filename, is_spherical, is_periodic, cell, translational_vector) AO_S += compute_overlaps(shell_2,shell_2p, nprocs) AO_St += compute_overlaps(shell_1,shell_2p, nprocs) + print('Done with computing atomic orbital overlaps. Elapsed time:', time.time()-t1) t1 = time.time() print('Turning the MATRIX to numpy array...') @@ -591,6 +606,7 @@ def run_cp2k_libint_step2(params): t1 = time.time() new_indices = CP2K_methods.resort_molog_eigenvectors(l_vals) eigenvectors_2 = [] + for j in range(len(eig_vect_2)): # the new and sorted eigenvector eigenvector_2 = eig_vect_2[j] @@ -611,12 +627,13 @@ def run_cp2k_libint_step2(params): ## t1 = time.time() print('Computing and saving molecular orbital overlaps...') + if isUKS: S_alpha = np.linalg.multi_dot([alpha_eigenvectors_2, AO_S, alpha_eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] - St_alpha = np.linalg.multi_dot([alpha_eigenvectors_1, AO_S, alpha_eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] + St_alpha = np.linalg.multi_dot([alpha_eigenvectors_1, AO_St, alpha_eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] S_beta = np.linalg.multi_dot([beta_eigenvectors_2, AO_S, beta_eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] - St_beta = np.linalg.multi_dot([beta_eigenvectors_1, AO_S, beta_eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] + St_beta = np.linalg.multi_dot([beta_eigenvectors_1, AO_St, beta_eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] S_step = data_conv.form_block_matrix(S_alpha,zero_mat,zero_mat.T,S_beta) St_step = data_conv.form_block_matrix(St_alpha,zero_mat,zero_mat.T,St_beta) @@ -630,7 +647,7 @@ def run_cp2k_libint_step2(params): E_step_sparse = scipy.sparse.csc_matrix(E_step) else: S = np.linalg.multi_dot([eigenvectors_2, AO_S, eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] - St = np.linalg.multi_dot([eigenvectors_1, AO_S, eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] + St = np.linalg.multi_dot([eigenvectors_1, AO_St, eigenvectors_2.T])[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital] S_step = data_conv.form_block_matrix(S,zero_mat,zero_mat,S) St_step = data_conv.form_block_matrix(St,zero_mat,zero_mat,St) S_step_sparse = scipy.sparse.csc_matrix(S_step) @@ -639,21 +656,27 @@ def run_cp2k_libint_step2(params): zero_mat,zero_mat,\ np.diag(energies_2)[lowest_orbital-1:highest_orbital,lowest_orbital-1:highest_orbital]) E_step_sparse = scipy.sparse.csc_matrix(E_step) + scipy.sparse.save_npz(params['res_dir']+F'/S_ks_{step}.npz', S_step_sparse) scipy.sparse.save_npz(params['res_dir']+F'/St_ks_{step-1}.npz', St_step_sparse) scipy.sparse.save_npz(params['res_dir']+F'/E_ks_{step}.npz', E_step_sparse) + print('Done with computing molecular orbital overlaps. Elapsed time:', time.time()-t1) + shell_1 = shell_2 energies_1 = energies_2 eigenvectors_1 = eigenvectors_2 + if params['remove_molden']: os.system(F'rm {molden_filename}') + if isxTB: print('Removing unnecessary wfn files...') os.system(F'rm OT_{step-1}-RESTART*') os.system(F'rm Diag_{step-1}-RESTART*') else: os.system(F'rm Diag_libra-{step-1}-RESTART*') + if params['cube_visualization']: print('Plotting cube files using VMD...') num_orbitals = E_step.shape[0] @@ -666,24 +689,28 @@ def run_cp2k_libint_step2(params): # Beta orbitals cube_file_name = F'Diag_{step-1}-WFN_{str(state_to_plot).zfill(5)}_2-1_0.cube' cube_file_methods.plot_cube_v2(params, cube_file_name, phase_factors_beta[c_plot,0]) + if params['plot_phase_corrected']: - if St_step[state_index, state_index]>0: + if St_step[state_index, state_index] > 0: f_alpha = 1 else: f_alpha = -1 + if St_step[state_index+num_orbitals, state_index+num_orbitals]>0: f_beta = 1 else: f_beta = -1 phase_factors_alpha[c_plot,0] *= f_alpha phase_factors_alpha[c_plot,0] *= f_beta + else: # Only alpha orbitals state_index = state_to_plot-params['lowest_orbital'] cube_file_name = F'Diag_{step-1}-WFN_{str(state_to_plot).zfill(5)}_1-1_0.cube' cube_file_methods.plot_cube_v2(params, cube_file_name, phase_factors_alpha[c_plot,0]) + if params['plot_phase_corrected']: - if St_step[state_index, state_index]>0: + if St_step[state_index, state_index] > 0: f_alpha = 1 else: f_alpha = -1 diff --git a/src/libra_py/workflows/nbra/step3.py b/src/libra_py/workflows/nbra/step3.py index 0bfb4d08a..0e706116e 100755 --- a/src/libra_py/workflows/nbra/step3.py +++ b/src/libra_py/workflows/nbra/step3.py @@ -1981,7 +1981,8 @@ def run_step3_ks_nacs_libint(params): St_ks = np.array( sp.load_npz(F'{res_dir_2}/St_ks_orthonormalized_{step}.npz').todense().real ) St_ks_cmatrix = data_conv.nparray2CMATRIX(St_ks) St_ks_cmatrices.append(St_ks_cmatrix) - E_ks = np.array( sp.load_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_re.npz').todense().real ) + ###E_ks = np.array( sp.load_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_re.npz').todense().real ) + E_ks = np.array( sp.load_npz(F'{res_dir_2}/E_ks_{step}.npz').todense().real ) E_ks_cmatrix = data_conv.nparray2CMATRIX(E_ks) E_ks_cmatrices.append(E_ks_cmatrix) # Now apply the state-reordering @@ -1992,7 +1993,8 @@ def run_step3_ks_nacs_libint(params): St_ks_sparse = data_conv.MATRIX2scipynpz( St_ks_cmatrices[step-start_time].real() ) sp.save_npz(F'{res_dir_2}/St_ks_orthonormalized_{step}.npz', St_ks_sparse) E_ks_sparse = data_conv.MATRIX2scipynpz( E_ks_cmatrices[step-start_time].real() ) - sp.save_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_re.npz', E_ks_sparse) + ###sp.save_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_re.npz', E_ks_sparse) + sp.save_npz(F'{res_dir_2}/E_ks_{step}.npz', E_ks_sparse) print('Done with state-reordering of the KS orbitals. Elapsed time:',time.time()-t2) # Applying phase correction @@ -2012,20 +2014,23 @@ def run_step3_ks_nacs_libint(params): St_step_phase_corrected, cum_phase_aa, cum_phase_bb = \ apply_phase_correction_scipy(St_step, step, cum_phase_aa, cum_phase_bb, two_spinor_format=True) Hvib_ks = 0.5/dt * (St_step_phase_corrected.todense().T - St_step_phase_corrected.todense()) - sp.save_npz(F'{res_dir_2}/St_ks_{step-start_time}_re.npz', St_step_phase_corrected ) - sp.save_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_im.npz', sp.csc_matrix( Hvib_ks )) + ###sp.save_npz(F'{res_dir_2}/St_ks_{step-start_time}_re.npz', St_step_phase_corrected ) + ###sp.save_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_im.npz', sp.csc_matrix( Hvib_ks )) + sp.save_npz(F'{res_dir_2}/St_ks_{step}.npz', St_step_phase_corrected ) + sp.save_npz(F'{res_dir_2}/Hvib_ks_{step}_im.npz', sp.csc_matrix( Hvib_ks )) os.system(F'rm {res_dir_2}/St_ks_orthonormalized_{step}.npz') else: for step in range(start_time, finish_time): St_step = sp.load_npz(F'{res_dir_2}/St_ks_orthonormalized_{step}.npz') Hvib_ks = 0.5/dt * (St_step.todense().T - St_step.todense()) - sp.save_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_im.npz', sp.csc_matrix( Hvib_ks )) - os.system(F'mv {res_dir_2}/St_ks_orthonormalized_{step}.npz {res_dir_2}/St_ks_{step-start_time}_re.npz') + ###sp.save_npz(F'{res_dir_2}/Hvib_ks_{step-start_time}_im.npz', sp.csc_matrix( Hvib_ks )) + sp.save_npz(F'{res_dir_2}/Hvib_ks_{step}_im.npz', sp.csc_matrix( Hvib_ks )) + ###os.system(F'mv {res_dir_2}/St_ks_orthonormalized_{step}.npz {res_dir_2}/St_ks_{step-start_time}_re.npz') + os.system(F'mv {res_dir_2}/St_ks_orthonormalized_{step}.npz {res_dir_2}/St_ks_{step}.npz') print('Done with phase correction. Elapsed time:', time.time()-t2) - def orthonormalize_ks_overlaps(step, params): """ This fuction is an auiliary function used in run_step3_ks_nacs_libint function @@ -2055,7 +2060,9 @@ def orthonormalize_ks_overlaps(step, params): S_step_sparse = sp.csc_matrix(S_step) sp.save_npz(F'{params["path_to_save_ks_Hvibs"]}/St_ks_orthonormalized_{step}.npz', St_step_sparse) - sp.save_npz(F'{params["path_to_save_ks_Hvibs"]}/Hvib_ks_{step-start_time}_re.npz', E_step) + sp.save_npz(F'{params["path_to_save_ks_Hvibs"]}/S_ks_{step}.npz', S_step_sparse) + ###sp.save_npz(F'{params["path_to_save_ks_Hvibs"]}/Hvib_ks_{step-start_time}_re.npz', E_step) + sp.save_npz(F'{params["path_to_save_ks_Hvibs"]}/E_ks_{step}.npz', E_step) print('Done with step', step,'. Elapsed time:', time.time()-t2) @@ -2147,43 +2154,58 @@ def run_step3_sd_nacs_libint(params): if params['is_many_body']: params['isnap'] = params['start_time'] params['fsnap'] = params['finish_time'] + #params['orbital_indices'] #ks_orbital_indicies = params['orbital_indices'] # The KS HOMO index specified by the user ks_homo_index = params['homo_index'] + # Generate the excitation data res = step3_many_body.get_step2_mb_sp_properties( params ) + # The unique SDs in the TD-DFT results in all logfiles sd_unique_basis = res[0] + # The ci_basis_states ci_basis_states = res[1] + # The CI coefficients ci_coefficients = res[2] + # The TD-DFT excitation energies ci_energies = res[3] + # And their spin-components (alpha and beta) spin_components = res[4] + # Now we need to generate the active space and KS HOMO index in that # active space for use in reading the npz files # The target KS states in TD-DFT sd_fstates = [] + # The initial KS states in TD-DFT sd_tstates = [] for i in range(len(sd_unique_basis)): sd_fstates.append(sd_unique_basis[i][0][0]) sd_tstates.append(sd_unique_basis[i][0][1]) + # The min_band present in the excitation min_band = min(sd_fstates) + # The max_band present in the excitation max_band = max(sd_tstates) + # number of occupied and unoccupied KS states num_occ = ks_homo_index-min_band+1 num_unocc = max_band-ks_homo_index+1 + # The new KS active space and HOMO index ks_active_space, ks_homo_index_1 = make_active_space(num_occ, num_unocc, params['data_dim'], params['npz_file_ks_homo_index']) params['active_space'] = ks_active_space + # The KS orbital indices ks_orbital_indicies = range(min_band, max_band+1) + # This parameter will be used in the reindexing of the SDs params['orbital_indices'] = ks_orbital_indicies #ks_homo_index = params['homo_index'] @@ -2191,11 +2213,13 @@ def run_step3_sd_nacs_libint(params): # Create single-particle bases based on the # the number of occupied and unoccupied KS states sd_unique_basis = [] + # In this case the min_band will be 1 and the # max_band is the data_dim/2 min_band = 1 max_band = int(data_dim/2) ks_orbital_indicies = range(min_band, max_band+1) + # Create the new active space and generate the KS HOMO index in that active space ks_active_space, ks_homo_index = make_active_space(params['num_occ_states'], params['num_unocc_states'], params['data_dim'], params['npz_file_ks_homo_index']) @@ -2207,17 +2231,20 @@ def run_step3_sd_nacs_libint(params): max_band = ks_homo_index+params['num_unocc_states'] for occ_state in range(ks_homo_index-params['num_occ_states']+1, ks_homo_index+1): for unocc_state in range(ks_homo_index+1,max_band+1): + # Alpha spin excitations sd_tmp = [[occ_state, unocc_state], 'alp'] if sd_tmp not in sd_unique_basis: sd_unique_basis.append(sd_tmp) - # If unrestricted spin calculations is requested the - # add the beta spin excitations as well + + # If unrestricted spin calculations is requested the + # add the beta spin excitations as well if params['isUKS']==1: sd_tmp = [[occ_state, unocc_state], 'bet'] #sd_tmp = [[state, ks_homo_index], 'bet'] if sd_tmp not in sd_unique_basis: sd_unique_basis.append(sd_tmp) + # I will keep this part for now as previously used and only commented the hole-only part ## Add hole-only SDs min_band = ks_homo_index-params['num_occ_states']+1 @@ -2230,8 +2257,10 @@ def run_step3_sd_nacs_libint(params): # sd_tmp = [[state, ks_homo_index+1], 'bet'] # if sd_tmp not in sd_unique_basis: # sd_unique_basis.append(sd_tmp) + # The params is updated with sd_unique_basis for reindexing params['sd_unique_basis'] = sd_unique_basis + # Reindex the SD basis sd_states_reindexed = step2_many_body.reindex_cp2k_sd_states( ks_homo_index, ks_orbital_indicies, sd_unique_basis, sd_format=2 ) @@ -2242,6 +2271,7 @@ def run_step3_sd_nacs_libint(params): print('ks_orbital_indicies', ks_orbital_indicies) params['isnap'] = start_time params['fsnap'] = finish_time + # The flag for phase-correction algorithm apply_phase_correction = params['apply_phase_correction'] @@ -2250,17 +2280,20 @@ def run_step3_sd_nacs_libint(params): E_sds = [] sd_states_reindexed_sorted = [] sd_states_unique_sorted = [] + for step in range(start_time, finish_time): E_sd_step, sd_states_unique_sorted_step, sd_states_reindexed_sorted_step, reindex_nsteps = sort_unique_SD_basis_scipy( step, sd_unique_basis, sd_states_reindexed, params ) sd_states_reindexed_sorted.append(sd_states_reindexed_sorted_step) sd_states_unique_sorted.append(sd_states_unique_sorted_step) E_sds.append(np.array(E_sd_step.todense().real)) - if step>start_time: + + if step > start_time: E_midpoint = 0.5*(E_sd_step+E_sd_step_plus) sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/Hvib_sd_{step-1}_re.npz',E_midpoint) E_sd_step_plus = E_sd_step print('Done with sorting and computing the SDs energies. Elapsed time:',time.time()-t2) + # Update the params with sd_states_reindexed_sorted params.update({'sd_states_reindexed_sorted':sd_states_reindexed_sorted}) if params['is_many_body']: @@ -2273,10 +2306,12 @@ def run_step3_sd_nacs_libint(params): else: # The old way SD2CI = step3_many_body.make_T_matrices( ci_coefficients, ci_basis_states, spin_components, sd_states_unique_sorted, params ) + # These parameters are in CMATRIX format. We need # to convert them into numpy arrays for i in range(len(SD2CI)): SD2CI[i] = data_conv.MATRIX2nparray(SD2CI[i].real()) + # We create the var_pool whether we use the multiprocessing or not var_pool = [] for step in range( finish_time - start_time -1 ): @@ -2297,6 +2332,7 @@ def run_step3_sd_nacs_libint(params): for variable in var_pool: St_sd = compute_sd_overlaps_in_parallel(variable[0],variable[1]) St_sds.append(St_sd) + print('Done with computing the SD overlaps. Elapsed time:', time.time()-t2) if params['do_state_reordering']==2 or params['do_state_reordering']==1: t2 = time.time() @@ -2308,9 +2344,11 @@ def run_step3_sd_nacs_libint(params): St_sds_cmatrix.append(St_sd_cmatrix) E_sd_cmatrix = data_conv.nparray2CMATRIX(E_sds[i]) E_sds_cmatrix.append(E_sd_cmatrix) + print('Applying state-reordering to SDs overlaps...\n') apply_state_reordering_general(St_sds_cmatrix, E_sds_cmatrix, params) print('Done with state-reordering of SDs. Elapsed time:', time.time()-t2) + if params['is_many_body']: St_cis = [] if params['do_state_reordering']==2 or params['do_state_reordering']==1: @@ -2323,12 +2361,14 @@ def run_step3_sd_nacs_libint(params): # Compute the St_ci St_ci = np.linalg.multi_dot([sd2ci.T, St_sds[i].todense().real, sd2ci]) St_cis.append(sp.csc_matrix(St_ci)) + # Now we need to apply state-reordering to St_cis # Again the same procedure is needed to convert between data types # Set up the counter c = 0 E_cis_cmatrix = [] St_cis_cmatrix = [] + for step in range(start_time, finish_time-1): # Adding 0.0 to the total energy E_ci_step = np.concatenate( (np.array([0.0]), np.array(ci_energies[step-start_time])) ) @@ -2339,10 +2379,12 @@ def run_step3_sd_nacs_libint(params): sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/Hvib_ci_{step-1}_re.npz',sp.csc_matrix(E_midpoint)) E_ci_step_plus = E_ci_step c += 1 + #if params['do_state_reordering']==2 or params['do_state_reordering']==1: print('Applying state-reordering to St_ci...\n') apply_state_reordering_general(St_cis_cmatrix, E_cis_cmatrix, params) print('Done with state-reordering to St_cis. Elapsed time:', time.time()-t2) + # Now convert back to scipy npz for i in range(len(St_cis_cmatrix)): St_cis[i] = data_conv.MATRIX2scipynpz( St_cis_cmatrix[i].real() ) @@ -2356,6 +2398,7 @@ def run_step3_sd_nacs_libint(params): # Now, we move to applying phase-corrections if params['apply_phase_correction']: + t3 = time.time() for step in range(len(St_sds)): # Initialize the nstates and cum_phase_aa and cum_phase_bb @@ -2364,12 +2407,15 @@ def run_step3_sd_nacs_libint(params): nstates = int(St_sds[0].shape[0]) #/2) cum_phase_aa = np.ones((nstates,1)) cum_phase_bb = np.ones((nstates,1)) + St_step_phase_corrected, cum_phase_aa, cum_phase_bb = \ apply_phase_correction_scipy(St_sds[step].real, step, cum_phase_aa, cum_phase_bb, two_spinor_format=False) + Hvib_sd = 0.5/dt * (St_step_phase_corrected.todense().T - St_step_phase_corrected.todense()) sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/Hvib_sd_{step+start_time}_im.npz', sp.csc_matrix( Hvib_sd )) sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/St_sd_{step+start_time}_re.npz', St_step_phase_corrected ) print('Done with applying phase-correction to St_sd matrices. Elpased time:', time.time()-t3) + if params['is_many_body']: # Set up the counter c = 0 @@ -2380,8 +2426,10 @@ def run_step3_sd_nacs_libint(params): nstates = int(St_cis[0].shape[0]) cum_phase_aa = np.ones((nstates,1)) cum_phase_bb = np.ones((nstates,1)) + St_ci_step_phase_corrected, cum_phase_aa, cum_phase_bb = \ apply_phase_correction_scipy(St_cis[step].real, step, cum_phase_aa, cum_phase_bb, two_spinor_format=False) + Hvib_ci = 0.5/dt * (St_ci_step_phase_corrected.todense().T - St_ci_step_phase_corrected.todense()) sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/Hvib_ci_{step+start_time}_im.npz', sp.csc_matrix( Hvib_ci )) sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/St_ci_{step+start_time}_re.npz', sp.csc_matrix(St_ci_step_phase_corrected)) @@ -2393,6 +2441,7 @@ def run_step3_sd_nacs_libint(params): Hvib_sd = 0.5/dt * (St_sds[step].todense().real.T - St_sds[step].todense().real) sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/Hvib_sd_{step+start_time}_im.npz', sp.csc_matrix( Hvib_sd )) sp.save_npz(F'{params["path_to_save_sd_Hvibs"]}/St_sd_{step+start_time}_re.npz', St_sds[step] ) + if params['is_many_body']: for step in range( finish_time - start_time -1 ): St_ci_step = St_cis[step].todense().real @@ -2403,17 +2452,19 @@ def run_step3_sd_nacs_libint(params): def compute_sd_overlaps_in_parallel( step, params ): - """ This function is used as an auxilary function that computes the SDs overlaps. It is used in the run_step3_sd_nacs_libint function. """ + start_time = params['start_time'] res_dir_1 = params['path_to_npz_files'] active_space = params['active_space'] sd_states_reindexed_sorted = params['sd_states_reindexed_sorted'] + t2 = time.time() print('Computing the SD overlaps for step',step) + # The same procedure as above st_ks = np.array( sp.load_npz(F'{res_dir_1}/St_ks_{step+start_time}.npz').todense() [active_space,:][:,active_space] ).real @@ -2421,9 +2472,11 @@ def compute_sd_overlaps_in_parallel( step, params ): [active_space,:][:,active_space] ).real s_ks_2 = np.array( sp.load_npz(F'{res_dir_1}/S_ks_{step+start_time}.npz').todense() [active_space,:][:,active_space] ).real + st_ks = data_conv.nparray2MATRIX(st_ks) s_ks_1 = data_conv.nparray2MATRIX(s_ks_1) s_ks_2 = data_conv.nparray2MATRIX(s_ks_2) + # Computing the overlaps for SDs t2 = time.time() s_sd_1 = mapping.ovlp_mat_arb(sd_states_reindexed_sorted[step], @@ -2432,10 +2485,13 @@ def compute_sd_overlaps_in_parallel( step, params ): sd_states_reindexed_sorted[step+1], s_ks_2,use_minimal=False, use_mo_approach=True).real() st_sd = mapping.ovlp_mat_arb(sd_states_reindexed_sorted[step], sd_states_reindexed_sorted[step+1], st_ks, use_minimal=False, use_mo_approach=True).real() + s_sd_1 = data_conv.MATRIX2nparray(s_sd_1) s_sd_2 = data_conv.MATRIX2nparray(s_sd_2) st_sd = data_conv.MATRIX2nparray(st_sd) + + if params['apply_orthonormalization']: print('Applying orthonormalization for SDs for step', step) st_sd, s_sd = apply_orthonormalization_scipy(s_sd_1.real, s_sd_2.real, st_sd.real)