From abc4b5d7ca2a210712c38ad625415efd5829aaf9 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 4 Dec 2023 10:05:55 -0500 Subject: [PATCH 01/47] Added a revised update_proj_adi function and used it in the main dynamics routine, also use the new active state update function there, added a couple of new experimental electronic integration schemes --- src/dyn/Dynamics.cpp | 124 +++++++++++++++++++++++++++++++++++-------- 1 file changed, 103 insertions(+), 21 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 6267284e2..e9aced393 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -570,6 +570,38 @@ void remove_thermal_correction(dyn_variables& dyn_var, nHamiltonian& ham, dyn_co } +/* +void update_proj_adi(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonian* Ham_prev, dyn_control_params& prms){ + + Just re-compute the proj_adi matrices + + + + //======= Parameters of the dyn variables ========== + int ntraj = dyn_var.ntraj; + + for(int itraj=0; itraj=100 && method <200){ traj1 = 0; } + + nHamiltonian* ham = Ham->children[traj1]; + nHamiltonian* ham_prev = Ham_prev->children[traj1]; + + //================= Basis re-expansion =================== + CMATRIX P(ham->nadi, ham->nadi); + CMATRIX T(*dyn_var.proj_adi[itraj]); T.load_identity(); + CMATRIX T_new(ham->nadi, ham->nadi); + + P = ham->get_time_overlap_adi(); // U_old.H() * U; + + // More consistent with the new derivations: + FullPivLU_inverse(P, T_new); + T_new = orthogonalized_T( T_new ); + + *dyn_var.proj_adi[itraj] = T_new; + }//for ntraj + +}// reproject_basis +*/ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonian* Ham_prev, dyn_control_params& prms){ @@ -591,6 +623,8 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia if(prms.rep_tdse==0 || prms.rep_tdse==2 ){ nst = ndia; } else if(prms.rep_tdse==1 || prms.rep_tdse==3 ){ nst = nadi; } + int ampl_transformation_method = prms.ampl_transformation_method; + CMATRIX C(nst, 1); CMATRIX vRHO(nst*nst, 1); // vectorized DM @@ -611,7 +645,7 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia nHamiltonian* ham = Ham->children[traj1]; nHamiltonian* ham_prev = Ham_prev->children[traj1]; - +/* //================= Basis re-expansion =================== CMATRIX P(ham->nadi, ham->nadi); @@ -629,7 +663,14 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia FullPivLU_inverse(P, T_new); T_new = orthogonalized_T( T_new ); - if(prms.assume_always_consistent){ T_new.identity(); } + *dyn_var.proj_adi[itraj] = T_new; + +*/ + CMATRIX T_new(*dyn_var.proj_adi[itraj]); + if(prms.assume_always_consistent){ T_new.identity(); } + + +// if ampl_transformation_method if(rep==0){ // diabatic @@ -678,6 +719,7 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia 3 - 1-point, Hvib integration, with exp_ 4 - 2-points, Hvib integration, with exp_ 5 - 2-points, Hvib, integration with the second-point correction of Hvib, with exp_ + 6 - 2-points, Hvib, no reordering 10 - same as 0, but with rotations 11 - same as 1, but with rotations @@ -731,6 +773,7 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia A = exp_(Hvib, complex(0.0, -0.5*dt) ); C = T_new * A * C; +// C = A * C; }// method == 2 && 102 @@ -785,7 +828,49 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia A = exp_(Hvib, complex(0.0, -0.5*dt) ); C = A * C; //T_new.identity(); - }// method == 4 && 104 + }// method == 6 && 106 + + else if(method==7 || method==107){ + // 2-point Hvib with vibronic Hamiltonian and correction - same as 5, but different order + // Propagate old coefficients using old Hamiltonian, only half-way + Hvib_old = ham_prev->get_hvib_adi(); // T is the identity matrix + if(is_ssy){ SSY_correction(Hvib_old, dyn_var, ham_prev, itraj); } + A = exp_(Hvib_old, complex(0.0, -0.5*dt) ); + C = A * C; + + // Reproject the coefficients + C = T_new.H() * C; + + // Propagate the new coefficients using the transformed new Hamiltonian, half-way + Hvib = ham->get_hvib_adi(); + Hvib = T_new.H() * Hvib * T_new; + if(is_ssy){ SSY_correction(Hvib, dyn_var, ham, itraj); } + + A = exp_(Hvib, complex(0.0, -0.5*dt) ); + C = A * C; + + }// method == 7 && 107 + + else if(method==8 || method==108){ + // new LD: 2-point Ham - same as 2 but different order + + // Propagate with the old Ham + Hvib_old = ham_prev->get_ham_adi(); // T is the identity matrix + if(is_ssy){ SSY_correction(Hvib_old, dyn_var, ham_prev, itraj); } + A = exp_(Hvib, complex(0.0, -0.5*dt) ); + C = A * C; + + // Reorder: + C = T_new.H() * C; + + // Propagate with the new Ham + Hvib = ham->get_ham_adi(); + if(is_ssy){ SSY_correction(Hvib, dyn_var, ham, itraj); } + Hvib = T_new.H() * Hvib * T_new; + A = exp_(Hvib, complex(0.0, -0.5*dt) ); + C = A * C; + + }// method == 8 && 108 @@ -972,7 +1057,7 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia - *dyn_var.proj_adi[itraj] = T_new; +// *dyn_var.proj_adi[itraj] = T_new; // Insert the propagated result back for(int st=0; st Date: Mon, 4 Dec 2023 10:09:34 -0500 Subject: [PATCH 02/47] Added new update_active_states function, redefined the meaning of some input variables to it, made more flexible, revised the internals, added some math documentation, also added some theory documentation to the orthgonalized_T function --- src/dyn/dyn_variables_electronic.cpp | 125 ++++++++++++++++++++------- 1 file changed, 96 insertions(+), 29 deletions(-) diff --git a/src/dyn/dyn_variables_electronic.cpp b/src/dyn/dyn_variables_electronic.cpp index 77ee80c09..7e0aa366b 100644 --- a/src/dyn/dyn_variables_electronic.cpp +++ b/src/dyn/dyn_variables_electronic.cpp @@ -744,52 +744,119 @@ void dyn_variables::init_electronic_dyn_var(bp::dict _params, Random& rnd){ //vector update_active_states(vector& act_states, vector& T){ -void dyn_variables:: update_active_states(){ +void dyn_variables:: update_active_states(int direction, int property){ /* + |psi_adi_tilde> C_adi_tilde = |psi_adi> C_adi + + |psi_adi_tilde> = |psi_adi> T, so: + + T C_adi_tilde = C_adi and C_adi_tilde = T^+ C_adi + + direction = 1 - forward: C_adi_tilde = T^+ C_adi + -1 - backward C_adi = T C_adi_tilde + + property = 0 - active state only + = 1 - adiabatic amplitudes only + = 2 - both active states and adiabatic amplitudes + act_states[itraj] - index of the active adiabatic state for the trajectory `itraj` in terms of the energy-ordered (instantaneous) eign-states */ - //cout<<"========= in update_active_states=============\n"; - //int ntraj = act_states.size(); - //int nst = T[0]->n_cols; - //vector res(ntraj, 0); - + int i; int nst = nadi; - - CMATRIX t2(nst,nst); - CMATRIX t2_half(nst, nst); - CMATRIX t2_i_half(nst, nst); CMATRIX coeff(nst, 1); complex max_val; for(int itraj=0; itraj "<(1.0, 0.0)); + + if(direction==1){ coeff = proj_adi[itraj]->H() * coeff; } + else if(direction==-1){ coeff = (*proj_adi[itraj]) * coeff; } + + coeff.max_col_elt(0, max_val, act_states[itraj]); + } + + //============ Amplitudes of states =================== + if(property==1 || property==2){ + coeff = ampl_adi->col(itraj); + + if(direction==1){ coeff = proj_adi[itraj]->H() * coeff; } + else if(direction==-1){ coeff = (*proj_adi[itraj]) * coeff; } + + for(i=0; iset(i, itraj, coeff.get(i, 0) ); } + } }// for itraj +} + +void dyn_variables::update_active_states(){ +/** + for backward-compatability +*/ + update_active_states(1, 0); +} - //cout<<"============ done with the update ================\n"; -// return res; + + + + +/* +void dyn_variables::update_ampl_reproject(int option){ + + changes all the adiabatic properties of the Hamiltonian by the matrix T + + option = 0 - use matrix U = T + = 1 - use matrix U = T.H() + + + + CMATRIX U(*T); // option == 0 + if (option==1){ U = CMATRIX(T->H()); } + else if(option==2 || option==4){ + } +*/ + + + + + CMATRIX orthogonalized_T(CMATRIX& T){ -/* +/** + Return a unitary equivalent of the matrix T + + |psi_adi_tilde(t+dt)> = |psi_adi(t+dt)> T(t+dt) + + = T(t+dt) + (should be close to I) (call it P (t, t+dt) ) + + So: + + I = P(t, t+dt) * T(t+dt) + + So, ideally T(t+dt) should be P^{-1}, but we also want it to be orthonormal + + T = P^{-1} * N + + We want: T^+ T = N^+ (P^{-1})^+ P^(-1) N = I + + N = (X^+ X)^-{1/2}, where X = P^(-1), indeed: + + N^+ (P^{-1})^+ P^(-1) N = ((X^+ X)^{-1/2})^+ (X^+ X) (X^+ X )^{-1/2} = + + = [((X^+ X)^{-1/2})^+ ] [ (X^+ X )^{1/2}] = (X X^+)^{-1/2} (X^+ X)^{1/2} = + + = (X^+)^{-1/2} X^{-1/2} X^{1/2} (X^+)^{1/2} = (X^+)^{-1/2} (X^+)^{1/2} = I + + The input to this function should be P^{-1} + */ int nst = T.n_cols; From eb80184bc743475fc520a4d7f9e696dd4c2b4709 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 4 Dec 2023 10:14:03 -0500 Subject: [PATCH 03/47] Ehrenfest forces are now computed without T matrix transformation, basically using the raw adiabatic amplitudes and matrix elements, but those are made consistent at the integration step --- src/dyn/Energy_and_Forces.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/dyn/Energy_and_Forces.cpp b/src/dyn/Energy_and_Forces.cpp index a798cadb0..2b256faa9 100755 --- a/src/dyn/Energy_and_Forces.cpp +++ b/src/dyn/Energy_and_Forces.cpp @@ -316,7 +316,8 @@ void update_forces(dyn_control_params& prms, dyn_variables& dyn_vars, nHamiltoni prms.electronic_integrator==12 ){ option = 1; } option = 0; - *dyn_vars.f = ham.Ehrenfest_forces_adi(*dyn_vars.ampl_adi, 1, option, dyn_vars.proj_adi).real(); + //*dyn_vars.f = ham.Ehrenfest_forces_adi(*dyn_vars.ampl_adi, 1, option, dyn_vars.proj_adi).real(); + *dyn_vars.f = ham.Ehrenfest_forces_adi(*dyn_vars.ampl_adi, 1, option).real(); /* CMATRIX& U = *ham.basis_transform; CMATRIX& C = *dyn_vars.ampl_adi; From 16a944cca4ca934dae8927017e7dea0f4bd6a2af Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 4 Dec 2023 10:15:33 -0500 Subject: [PATCH 04/47] Added the ampl_transformation_method parameter but it is not currently used --- src/dyn/dyn_control_params.cpp | 3 +++ src/dyn/dyn_control_params.h | 8 ++++++++ 2 files changed, 11 insertions(+) diff --git a/src/dyn/dyn_control_params.cpp b/src/dyn/dyn_control_params.cpp index fabd5cec9..6078f3e3e 100755 --- a/src/dyn/dyn_control_params.cpp +++ b/src/dyn/dyn_control_params.cpp @@ -104,6 +104,7 @@ dyn_control_params::dyn_control_params(){ dt = 41.0; num_electronic_substeps = 1; electronic_integrator = 0; + ampl_transformation_method = 1; assume_always_consistent = 0; thermally_corrected_nbra = 0; @@ -184,6 +185,7 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){ dt = x.dt; num_electronic_substeps = x.num_electronic_substeps; electronic_integrator = x.electronic_integrator; + ampl_transformation_method = x.ampl_transformation_method; assume_always_consistent = x. assume_always_consistent; decoherence_rates = new MATRIX(x.decoherence_rates->n_rows, x.decoherence_rates->n_cols); @@ -378,6 +380,7 @@ void dyn_control_params::set_parameters(bp::dict params){ else if(key=="dt") { dt = bp::extract(params.values()[i]); } else if(key=="num_electronic_substeps") { num_electronic_substeps = bp::extract(params.values()[i]); } else if(key=="electronic_integrator"){ electronic_integrator = bp::extract(params.values()[i]); } + else if(key=="ampl_transformation_method"){ ampl_transformation_method = bp::extract(params.values()[i]); } else if(key=="assume_always_consistent"){ assume_always_consistent = bp::extract(params.values()[i]); } else if(key=="thermally_corrected_nbra"){ thermally_corrected_nbra = bp::extract(params.values()[i]); } diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index 4d78363ea..9a57385ac 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -695,6 +695,14 @@ class dyn_control_params{ */ int electronic_integrator; + + /** + Whether transform the amplitudes by the T transformation matrix + 0 - do not transform by the T matrix (naive, but potentially correct approach) + 1 - do transform it (as in LD, but maybe not needed if we directly transform basis) + */ + int ampl_transformation_method; + /** If set to True (1), we will force the reprojection matrix T_new to be the identity matrix. This effectively From 9add4d0f2f581ef228c96e4ee98fd161ba25253c Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 4 Dec 2023 10:17:37 -0500 Subject: [PATCH 05/47] implementation of the new update_proj_adi function, also added the update_time_overlaps function, in case we only want update that variable, but nothing else --- src/dyn/dyn_ham.cpp | 54 +++++++++++++++++++++++++++++++++++++++++++++ src/dyn/dyn_ham.h | 5 +++++ 2 files changed, 59 insertions(+) diff --git a/src/dyn/dyn_ham.cpp b/src/dyn/dyn_ham.cpp index 96357a7d4..16cb972a3 100644 --- a/src/dyn/dyn_ham.cpp +++ b/src/dyn/dyn_ham.cpp @@ -117,6 +117,7 @@ void update_Hamiltonian_variables(dyn_control_params& prms, dyn_variables& dyn_v exit(0); } + //============================== Time-overlaps ====================== /// Don't update time-overlaps if(prms.time_overlap_method==0){ ;; } // maybe it is already pre-computed and stored @@ -131,6 +132,7 @@ void update_Hamiltonian_variables(dyn_control_params& prms, dyn_variables& dyn_v } }// 1 + }// update_type == 0 @@ -237,6 +239,58 @@ void update_Hamiltonian_variables(bp::dict prms, dyn_variables& dyn_var, } +void update_time_overlaps(dyn_control_params& prms, dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev){ + + int nadi = ham.nadi; + int ntraj = ham.children.size(); + + /// Don't update time-overlaps + if(prms.time_overlap_method==0){ ;; } // maybe it is already pre-computed and stored + + /// Compute the time-overlap directly, using previous MO vectors + else if(prms.time_overlap_method==1){ // Explicitly compute it: + + CMATRIX st(nadi, nadi); + for(int traj=0; trajget_basis_transform().H() * ham.children[traj]->get_basis_transform(); + ham.children[traj]->set_time_overlap_adi_by_val(st); + } + }// 1 + +} + +void update_proj_adi(dyn_control_params& prms, dyn_variables& dyn_var, nHamiltonian& Ham){ // nHamiltonian& Ham_prev){ +/** + Just re-compute the proj_adi matrices +*/ + + //======= Parameters of the dyn variables ========== + int ntraj = dyn_var.ntraj; + + for(int itraj=0; itraj=100 && method <200){ traj1 = 0; } + + nHamiltonian* ham = Ham.children[traj1]; + //nHamiltonian* ham_prev = Ham_prev.children[traj1]; + + //================= Basis re-expansion =================== + CMATRIX P(ham->nadi, ham->nadi); + CMATRIX T_new(*dyn_var.proj_adi[itraj]); + + P = ham->get_time_overlap_adi(); // (U_old.H() * U).H(); + + + // More consistent with the new derivations: + FullPivLU_inverse(P, T_new); + T_new = orthogonalized_T( T_new ); + + *dyn_var.proj_adi[itraj] = T_new; + }//for ntraj + +}// reproject_basis + + + }// namespace libdyn diff --git a/src/dyn/dyn_ham.h b/src/dyn/dyn_ham.h index 4f22ff598..a7a19be95 100644 --- a/src/dyn/dyn_ham.h +++ b/src/dyn/dyn_ham.h @@ -45,6 +45,11 @@ void update_Hamiltonian_variables(bp::dict prms, dyn_variables& dyn_var, bp::object py_funct, bp::object model_params, int update_type); +void update_time_overlaps(dyn_control_params& prms, dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev); + +void update_proj_adi(dyn_control_params& prms, dyn_variables& dyn_var, nHamiltonian& Ham); // nHamiltonian& Ham_prev); + + }// namespace libdyn }// liblibra From a211b854484cf9c4038db68513b6cd878b8c9fb5 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 4 Dec 2023 10:18:53 -0500 Subject: [PATCH 06/47] interface of the new update_active_states function --- src/dyn/dyn_variables.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/dyn/dyn_variables.h b/src/dyn/dyn_variables.h index dcb998852..0e55a1d54 100755 --- a/src/dyn/dyn_variables.h +++ b/src/dyn/dyn_variables.h @@ -494,6 +494,7 @@ class dyn_variables{ void update_density_matrix(dyn_control_params& dyn_params, bp::object compute_model, bp::dict model_params, int lvl); void update_density_matrix(bp::dict dyn_params, bp::object compute_model, bp::dict model_params, int lvl); + void update_active_states(int direction, int property); void update_active_states(); void init_amplitudes(bp::dict params, Random& rnd); From 3a01a37bb154b313fd786ac62d69f6f87a0fba53 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 4 Dec 2023 10:20:34 -0500 Subject: [PATCH 07/47] Exposed the new update_active_states function of the dyn_variables class to Python, also exposed the ampl_transformation_method member of the dyn_control_params class to Python --- src/dyn/libdyn.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index 1160f2572..d39718bd2 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -122,6 +122,7 @@ void export_dyn_control_params_objects(){ .def_readwrite("dt", &dyn_control_params::dt) .def_readwrite("num_electronic_substeps", &dyn_control_params::num_electronic_substeps) .def_readwrite("electronic_integrator", &dyn_control_params::electronic_integrator) + .def_readwrite("ampl_transformation_method", &dyn_control_params::ampl_transformation_method) .def_readwrite("assume_always_consistent", &dyn_control_params::assume_always_consistent) .def_readwrite("thermally_corrected_nbra", &dyn_control_params::thermally_corrected_nbra) .def_readwrite("total_energy", &dyn_control_params::total_energy) @@ -191,6 +192,9 @@ void export_dyn_variables_objects(){ vector (dyn_variables::*expt_compute_kinetic_energies_v2)(vector& which_dofs) = &dyn_variables::compute_kinetic_energies; + void (dyn_variables::*expt_update_active_states_v1)(int direction, int property) = &dyn_variables::update_active_states; + void (dyn_variables::*expt_update_active_states_v2)() = &dyn_variables::update_active_states; + class_("dyn_variables",init()) .def("__copy__", &generic__copy__) @@ -267,7 +271,8 @@ void export_dyn_variables_objects(){ .def("update_density_matrix", expt_update_density_matrix_v3) .def("update_density_matrix", expt_update_density_matrix_v4) - .def("update_active_states", &dyn_variables::update_active_states) + .def("update_active_states", expt_update_active_states_v1) + .def("update_active_states", expt_update_active_states_v2) .def("init_amplitudes", &dyn_variables::init_amplitudes) .def("init_density_matrix", &dyn_variables::init_density_matrix) From 206419a379d60541b4ea87227744b4bd69853256 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 4 Dec 2023 10:22:50 -0500 Subject: [PATCH 08/47] Implementation of the new transform_basis member of the nHamiltonian class, also the implementation of the changes of the meaning of arguments to transform_all function of the nHamiltonian class --- src/nhamiltonian/nHamiltonian.h | 2 ++ src/nhamiltonian/nHamiltonian_basic.cpp | 43 +++++++++++++++++++++---- 2 files changed, 39 insertions(+), 6 deletions(-) diff --git a/src/nhamiltonian/nHamiltonian.h b/src/nhamiltonian/nHamiltonian.h index 459dfac51..c7fb79999 100644 --- a/src/nhamiltonian/nHamiltonian.h +++ b/src/nhamiltonian/nHamiltonian.h @@ -422,6 +422,8 @@ class nHamiltonian{ void transform_all(CMATRIX* T, int option); void transform_all(vector& T, int option); + void transform_basis(CMATRIX* T, int option); + void transform_basis(vector& T, int option); ///< In nHamiltonian_compute_diabatic.cpp diff --git a/src/nhamiltonian/nHamiltonian_basic.cpp b/src/nhamiltonian/nHamiltonian_basic.cpp index de4aecbe8..0aa9fa804 100644 --- a/src/nhamiltonian/nHamiltonian_basic.cpp +++ b/src/nhamiltonian/nHamiltonian_basic.cpp @@ -2132,13 +2132,16 @@ void nHamiltonian::transform_all(CMATRIX* T, int option){ /** changes all the adiabatic properties of the Hamiltonian by the matrix T - option = 0 - use matrix U = T - = 1 - use matrix U = T.H() + option = 1 - use matrix U = T |psi_adi_tilde> = |psi_adi> T + = -1 - use matrix U = T.H() |psi_adi> = |psi_adi_tilde>T^+ */ - CMATRIX U(*T); // option == 0 - if (option==1){ U = CMATRIX(T->H()); } + CMATRIX U(nadi, nadi); + if(option==1){ U = *T; } + else if(option==-1){ U = T->H(); } + +/* else if(option==2 || option==4){ int nst = nadi; @@ -2161,7 +2164,7 @@ void nHamiltonian::transform_all(CMATRIX* T, int option){ FullPivLU_inverse(U, invU); U = invU; } - +*/ //FullPivLU_inverse(t, U); // U = t^{-1} @@ -2179,7 +2182,8 @@ void nHamiltonian::transform_all(CMATRIX* T, int option){ // Time-overlaps - *time_overlap_adi = U.H() * (*time_overlap_adi) * U; +// *time_overlap_adi = U.H() * (*time_overlap_adi) * U; + *time_overlap_adi = (*time_overlap_adi) * U; // Derivatives for(int i=0; i& T, int option){ } +void nHamiltonian::transform_basis(CMATRIX* T, int option){ +/** + changes the diabatic-to-adiabatic basis transformation matrix by a matrix T + + option = 1 - use matrix U = T |psi_adi_tilde> = |psi_adi> T + = -1 - use matrix U = T.H() |psi_adi> = |psi_adi_tilde>T^+ + +*/ + + CMATRIX U(nadi, nadi); + if(option==1){ U = *T; } + else if(option==-1){ U = T->H(); } + + // Basis transform + *basis_transform = (*basis_transform) * U; + +} // transform_basis + + +void nHamiltonian::transform_basis(vector& T, int option){ +/** + changes the basis of all the children Hamiltonians by the matrices in T +*/ + + int ntraj = children.size(); + for(int i=0; itransform_basis(T[i], option); } +}// transform_basis }// namespace libnhamiltonian }// liblibra From 3acb311478dc859c436ea4eefc7e9e0a82f4d01d Mon Sep 17 00:00:00 2001 From: MohammadShakiba Date: Mon, 4 Dec 2023 16:18:29 -0500 Subject: [PATCH 09/47] Fied eigenvalue sorting for KS equations, added molden file writing function --- src/libra_py/data_conv.py | 16 ++++++++++++ src/libra_py/data_read.py | 29 ++++++++++++++++++++++ src/libra_py/molden_methods.py | 35 ++++++++++++++++++++++++++- src/libra_py/packages/cp2k/methods.py | 2 +- 4 files changed, 80 insertions(+), 2 deletions(-) diff --git a/src/libra_py/data_conv.py b/src/libra_py/data_conv.py index a2ddc8901..95dc85da1 100755 --- a/src/libra_py/data_conv.py +++ b/src/libra_py/data_conv.py @@ -578,3 +578,19 @@ def vasp_to_xyz(filename): f.close() +def upper_vector_to_symmetric_nparray(upper_vector, upper_indices, mat_shape): + """ + This function gets the upper triangular part of a matrix as a vector and retuns a symmetric matrix + Args: + upper_vector (nparray): The upper triangular of a matrix + upper_indices (nparray): The indices of the upper part of the matrix + mat_shape (tuple): The shape of the original numpy array + Returns: + matrix (nparray): The symmetric matix built based on the upper triangular matrix + """ + matrix = np.zeros(mat_shape) + matrix[upper_indices] = upper_vector + matrix = matrix + matrix.T - np.diag(matrix.diagonal()) + return matrix + + diff --git a/src/libra_py/data_read.py b/src/libra_py/data_read.py index 5a36b4a9e..687566a18 100755 --- a/src/libra_py/data_read.py +++ b/src/libra_py/data_read.py @@ -434,4 +434,33 @@ def get_all_data(params): return Hadi, Hvib, nac, St, nstates + +def read_atom_coords(filename, step): + """ + This function reads the trajectory xyz file and returns the atoms and coordiantes of step^th geometry. + Args: + filename (string): The path to the trajectory xyz file. + step (integer): The step to be extracted. + Returns: + atoms (list): A list that contains the atoms names. + coords (nparray): A numpy array that contains the X,Y, and Z coordinates of each atom. + """ + f = open(filename, 'r') + lines = f.readlines() + f.close() + natoms = int(lines[0].split()[0]) + coords = [] + start = step*(natoms+2)+2 + end = start+natoms + atoms = [] + for i in range(start, end): + tmp = lines[i].split() + #print(lines[i]) + atoms.append(tmp[0]) + x = float(tmp[1]) + y = float(tmp[2]) + z = float(tmp[3]) + coords.append([x,y,z]) + return atoms, np.array(coords) + diff --git a/src/libra_py/molden_methods.py b/src/libra_py/molden_methods.py index 92995b043..6dd1655d8 100755 --- a/src/libra_py/molden_methods.py +++ b/src/libra_py/molden_methods.py @@ -26,7 +26,7 @@ elif sys.platform=="linux" or sys.platform=="linux2": from liblibra_core import * -from libra_py import data_outs +from libra_py import data_outs, data_read from libra_py import units import math @@ -402,3 +402,36 @@ def index_reorder(l_val): # The indeices return np.array(new_order)-1 +def write_molden_file(filename, sample_molden_file, path_to_trajectory, step): + """ + This function writes a molden file for a geometry that contains the basis data. It doesn't write the + eigenvectors since they are not needed anymore. + Args: + filename (strin): The molden file name to be written + sample_molden_file (strin): The path to the sample molden file that contains the basis data + path_to_trajectory (strin): The path to MD trajectory xyz file + step (integer): The MD time step + Returns: + None + """ + f = open(sample_molden_file) + lines = f.readlines() + f.close() + + for i in range(len(lines)): + if '[GTO]' in lines[i]: + start = i + elif '[5D7F]' in lines[i]: + end = i + + atoms, coords = data_read.read_atom_coords(path_to_trajectory, step) + coords = coords*units.Angst + f = open(filename,'w') + f.write('[Molden Format]\n [Atoms] AU\n') + for i in range(len(coords)): + f.write(F'{atoms[i]} {i+1} 0 {coords[i][0]} {coords[i][1]} {coords[i][2]}\n') + for i in range(start, end+1): + f.write(lines[i]) + f.close() + + diff --git a/src/libra_py/packages/cp2k/methods.py b/src/libra_py/packages/cp2k/methods.py index 65e285446..e1c59a68a 100755 --- a/src/libra_py/packages/cp2k/methods.py +++ b/src/libra_py/packages/cp2k/methods.py @@ -2020,7 +2020,7 @@ def compute_energies_coeffs(ks_mat, overlap): eigenvectors = eigenvectors[:,sorted_indices].T - return eigenvalues, eigenvectors + return eigenvalues[sorted_indices], eigenvectors def compute_density_matrix(eigenvectors, homo_index): """ From cd525631873e20a60e673176b7b8a1306ef47758 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Mon, 11 Dec 2023 12:15:56 -0500 Subject: [PATCH 10/47] Added new options for instantaneous decoherence methods: new approach, IDN, 3 and ID at failed transitions - IDF, option 4 --- src/dyn/Dynamics.cpp | 26 +++++- src/dyn/dyn_control_params.h | 8 +- src/dyn/dyn_decoherence.h | 3 + src/dyn/dyn_decoherence_methods.cpp | 134 +++++++++++++++++++++++++++- src/dyn/libdyn.cpp | 4 + 5 files changed, 168 insertions(+), 7 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index e9aced393..2c296b762 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -825,8 +825,11 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia Hvib += Hvib_old; - A = exp_(Hvib, complex(0.0, -0.5*dt) ); + A = exp_(Hvib_old, complex(0.0, -0.5*dt) ); C = A * C; + + //dyn_var.proj_adi[itraj]->load_identity(); + //T_new.identity(); }// method == 6 && 106 @@ -855,21 +858,24 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia // new LD: 2-point Ham - same as 2 but different order // Propagate with the old Ham - Hvib_old = ham_prev->get_ham_adi(); // T is the identity matrix + Hvib_old = ham_prev->get_hvib_adi(); // T is the identity matrix if(is_ssy){ SSY_correction(Hvib_old, dyn_var, ham_prev, itraj); } A = exp_(Hvib, complex(0.0, -0.5*dt) ); C = A * C; // Reorder: - C = T_new.H() * C; + //C = T_new.H() * C; // Propagate with the new Ham - Hvib = ham->get_ham_adi(); + Hvib = ham->get_hvib_adi(); if(is_ssy){ SSY_correction(Hvib, dyn_var, ham, itraj); } Hvib = T_new.H() * Hvib * T_new; A = exp_(Hvib, complex(0.0, -0.5*dt) ); C = A * C; + // Reorder back: + //C = T_new * C; + }// method == 8 && 108 @@ -1442,6 +1448,17 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, else if(prms.decoherence_algo==4){ ;; } // MF-SD of Bedard-Hearn, Larsen, Schwartz, nothing to do here //else if(prms.decoherence_algo==5 or prms.decoherence_algo==6){ xf_hop_reset(dyn_var, act_states, old_states); } // SHXF or MQCXF else if(prms.decoherence_algo==5){ xf_hop_reset(dyn_var, act_states, old_states); } // SHXF or MQCXF + else if(prms.decoherence_algo==6){ ;; } // MQCXF, nothing to do here? + else if(prms.decoherence_algo==7){ ;; } // DISH rev 2023, nothing to do here + + // Experimental: instantaneous decoherence in diabatic basis + else if(prms.decoherence_algo==8){ + if(prms.rep_tdse==1){ + instantaneous_decoherence_dia(*dyn_var.ampl_adi, ham, act_states, prop_states, old_states, + prms.instantaneous_decoherence_variant, prms.collapse_option); + } + else{ cout<<"ERROR: Instantaneous Decoherence requires rep_tdse = 1\nExiting now...\n"; exit(0); } + }// algo == 6 } // DISH - the old one @@ -1452,6 +1469,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, }// DISH + //====================== Momenta adjustment after successful/frustrated hops =================== // Velocity rescaling: however here we may be changing velocities handle_hops_nuclear(dyn_var, ham, act_states, old_states, prms); diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index 9a57385ac..b7ef90fe9 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -389,6 +389,7 @@ class dyn_control_params{ - 5: SHXF of Min - 6: MQCXF - 7: DISH, rev2023 + - 8: diabatic IDA, experimental */ double decoherence_algo; @@ -464,8 +465,13 @@ class dyn_control_params{ only used with decoherence_algo == 1 - 0: ID-S - - 1: ID-A [default] + - 1: ID-A [default] - if the proposed hop is not successful, we project back to the initial state + if the proposed hop is accepted - we project onto that state - 2: ID-C - consistent ID - an experimental algorithm + - 3: ID-A, new: if the proposed hop is not successful, we project out the proposed states + if the proposed hop is accepted - we project onto that state + - 4: ID-F, new: if the proposed hop is not successful, we project out the proposed states + but we don't do anything if the hop is successful */ int instantaneous_decoherence_variant; diff --git a/src/dyn/dyn_decoherence.h b/src/dyn/dyn_decoherence.h index b3d4f7b00..943706d78 100755 --- a/src/dyn/dyn_decoherence.h +++ b/src/dyn/dyn_decoherence.h @@ -54,6 +54,9 @@ void instantaneous_decoherence(CMATRIX& Coeff, vector& accepted_states, vector& proposed_states, vector& initial_states, int instantaneous_decoherence_variant, int collapse_option); +void instantaneous_decoherence_dia(CMATRIX& Coeff, nHamiltonian& ham, + vector& accepted_states, vector& proposed_states, vector& initial_states, + int instantaneous_decoherence_variant, int collapse_option); CMATRIX afssh_dzdt(CMATRIX& dz, CMATRIX& Hvib, CMATRIX& F, CMATRIX& C, double mass, int act_state); void integrate_afssh_moments(CMATRIX& dR, CMATRIX& dP, CMATRIX& Hvib, CMATRIX& F, CMATRIX& C, double mass, int act_state, double dt, int nsteps); diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index f24dd101f..892072272 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -353,6 +353,138 @@ void instantaneous_decoherence(CMATRIX& Coeff, } + // ID-S + if(instantaneous_decoherence_variant==0){ + + for(traj = 0; traj < ntraj; traj++){ + + if(accepted_states[traj] != initial_states[traj]){ + collapse(Coeff, traj, accepted_states[traj], collapse_option); + } + }// traj + + }// ID-S + + // ID-A + else if(instantaneous_decoherence_variant==1){ + for(traj = 0; traj < ntraj; traj++){ + if(proposed_states[traj] != initial_states[traj]){ + // Only apply ID-A, if the proposed states are different from the original ones + + if(accepted_states[traj] == proposed_states[traj]){ + // Proposed hop is successful - collapse onto newly accepted state + collapse(Coeff, traj, accepted_states[traj], collapse_option); + } + else{ + // Proposed hop is not successful - collapse onto the original state + collapse(Coeff, traj, initial_states[traj], collapse_option); + } + } + }// traj + }// ID-A + + // ID-C + else if(instantaneous_decoherence_variant==2){ + + for(traj = 0; traj < ntraj; traj++){ + collapse(Coeff, traj, accepted_states[traj], collapse_option); + }// traj + + }// ID-C + + // ID-A, new version = IDN + else if(instantaneous_decoherence_variant==3){ + for(traj = 0; traj < ntraj; traj++){ + if(proposed_states[traj] != initial_states[traj]){ + // Only apply ID-A, if the proposed states are different from the original ones + + if(accepted_states[traj] == proposed_states[traj]){ + // Proposed hop is successful - collapse onto newly accepted state + collapse(Coeff, traj, accepted_states[traj], collapse_option); + } + else{ + // Proposed hop is not successful - project out the proposed state + project_out(Coeff, traj, proposed_states[traj]); + } + } + }// traj + }// IDN + + // ID-F, ID on the failed transition, experimental + else if(instantaneous_decoherence_variant==4){ + for(traj = 0; traj < ntraj; traj++){ + if(proposed_states[traj] != initial_states[traj]){ + // Only apply ID-F, if the proposed states are different from the original ones + + if(accepted_states[traj] != proposed_states[traj]){ + // Proposed hop is rejected - project out the proposed states, but don't collapse + // onto the accepted hop + project_out(Coeff, traj, proposed_states[traj]); + } + } + }// traj + }// IDF + +} + + +void instantaneous_decoherence_dia(CMATRIX& Coeff, nHamiltonian& ham, + vector& accepted_states, vector& proposed_states, vector& initial_states, + int instantaneous_decoherence_variant, int collapse_option){ + +/** + This is an experimental ID function which does all the same as the normal one does, + but in the diabatic basis: + + Args: + Coeff - the adiabatic amplitudes + ham - nHamiltonian object + + Two options of the algorithm are available: + ID-S: wavefunction amplitudes are collapsed only during the successful hops (onto the new state) + ID-A: wavefunction amplitudes are collapsed at every attempted hop + to the new state, if successful + to the old state, if now + + There collapsing options are controlled by the parameter instantaneous_decoherence_variant + + 0 - ID-S + 1 - ID-A + 2 - ID-C - consistent ID - an experimental algo + + In the "consistent" version, we lift the condition that the accepted/proposed states must be different from + the starting state - this option addresses a philosophical question - what if we say that no hops + is equivalent to the hop onto the current state? why should the "hopping" into the original state + by treated differently from the "actual hopping" to another state? So, in this version we + collapse onto the accpeted states, no matter if they are the result of a successfull or frustrated hop. + + The mechanism of the collapse event itself is controlled by the collapse_option parameter +*/ + int traj; + int ntraj = Coeff.n_cols; + + if(accepted_states.size()!=ntraj){ + cout<<"ERROR in ids: the sizes of the input variables Coeff and accepted_states are inconsistent\n"; + cout<<"Coeff.num_of_cols = = "<& accepted_states, vector& proposed_states, vector& initial_states, + int instantaneous_decoherence_variant, int collapse_option) = &instantaneous_decoherence; + def("instantaneous_decoherence_dia", expt_instantaneous_decoherence_dia_v1); void (*expt_wp_reversal_events_v1) From 7d102c2af4851c043585062495b8b254adc10355 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Tue, 12 Dec 2023 10:15:15 -0500 Subject: [PATCH 11/47] Adjust the XF propagation according to the LD scheme --- src/dyn/Dynamics.cpp | 6 ++---- src/dyn/dyn_decoherence_methods.cpp | 13 +++++++++---- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 2c296b762..788dbce71 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1223,7 +1223,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, propagate_electronic(dyn_var, ham, ham_aux, prms); if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ - propagate_half_xf(dyn_var, ham, prms, 1); + propagate_half_xf(dyn_var, ham, prms, 0); // 1 } } @@ -1446,9 +1446,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, }// AFSSH else if(prms.decoherence_algo==3){ ;; } // BCSH of Linjun Wang, nothing to do here else if(prms.decoherence_algo==4){ ;; } // MF-SD of Bedard-Hearn, Larsen, Schwartz, nothing to do here - //else if(prms.decoherence_algo==5 or prms.decoherence_algo==6){ xf_hop_reset(dyn_var, act_states, old_states); } // SHXF or MQCXF - else if(prms.decoherence_algo==5){ xf_hop_reset(dyn_var, act_states, old_states); } // SHXF or MQCXF - else if(prms.decoherence_algo==6){ ;; } // MQCXF, nothing to do here? + else if(prms.decoherence_algo==5 or prms.decoherence_algo==6){ xf_hop_reset(dyn_var, act_states, old_states); } // SHXF or MQCXF else if(prms.decoherence_algo==7){ ;; } // DISH rev 2023, nothing to do here // Experimental: instantaneous decoherence in diabatic basis diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 892072272..632ba598e 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1549,7 +1549,8 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa int traj1 = itraj; if(method >=100 && method <200){ traj1 = 0; } nHamiltonian* ham = Ham.children[traj1]; - + +/* //================= Basis re-expansion =================== CMATRIX P(nadi, nadi); CMATRIX T(*dyn_var.proj_adi[itraj]); T.load_identity(); @@ -1562,8 +1563,12 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa //libmeigen::FullPivLU_inverse(P, T_new); //T_new = orthogonalized_T( T_new ); // - //if(prms.assume_always_consistent){ T_new.identity(); } - +*/ + + CMATRIX T(*dyn_var.proj_adi[itraj]); T.load_identity(); + CMATRIX T_new(*dyn_var.proj_adi[itraj]); + if(prms.assume_always_consistent){ T_new.identity(); } + // Generate the half-time exponential operator CMATRIX Hxf_old(nadi, nadi); CMATRIX Hxf(nadi, nadi); @@ -1584,7 +1589,7 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa C = D * C; } - *dyn_var.proj_adi[itraj] = T_new; +// *dyn_var.proj_adi[itraj] = T_new; // Insert the propagated result back for(int st=0; st Date: Tue, 12 Dec 2023 10:44:25 -0500 Subject: [PATCH 12/47] Remove redundant transformation of the XF Hamiltonian --- src/dyn/dyn_decoherence_methods.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 632ba598e..9dcb23b8d 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1575,7 +1575,7 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa CMATRIX D(nadi, nadi); /// this is \exp[-idt/4\hbar * ( T_new.H()*Hxf(t+dt)*T_new + Hxf(t) )] XF_correction(Hxf_old, dyn_var, C, prms.wp_width, T, itraj); - XF_correction(Hxf, dyn_var, C, prms.wp_width, T_new, itraj); + XF_correction(Hxf, dyn_var, C, prms.wp_width, T, itraj); Hxf = T_new.H() * Hxf * T_new; Hxf += Hxf_old; From a0c0c298dc5ff0d4a5d3775fb2dabe9460fc28ed Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Tue, 12 Dec 2023 18:20:42 -0500 Subject: [PATCH 13/47] Trim the code --- src/dyn/Dynamics.cpp | 4 ++-- src/dyn/dyn_decoherence.h | 2 +- src/dyn/dyn_decoherence_methods.cpp | 24 ++---------------------- 3 files changed, 5 insertions(+), 25 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 788dbce71..6fc13cd8a 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1217,13 +1217,13 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, // basis re-projection matrices for(i=0; iget_time_overlap_adi(); // U_old.H() * U; - // - //// More consistent with the new derivations: - //libmeigen::FullPivLU_inverse(P, T_new); - //T_new = orthogonalized_T( T_new ); - // -*/ - CMATRIX T(*dyn_var.proj_adi[itraj]); T.load_identity(); CMATRIX T_new(*dyn_var.proj_adi[itraj]); if(prms.assume_always_consistent){ T_new.identity(); } @@ -1582,12 +1567,7 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa D = libspecialfunctions::exp_(Hxf, complex(0.0, -0.25*dt) ); - if(do_rotation==1){ - C = T_new * D * T_new.H() * C; - } - else{ - C = D * C; - } + C = D * C; // *dyn_var.proj_adi[itraj] = T_new; From b47d81ff18b05f598a10a54a30c5d699d139f0c6 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Wed, 13 Dec 2023 07:57:51 -0500 Subject: [PATCH 14/47] Edit printout in MQCXF --- src/dyn/dyn_decoherence_methods.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 959092c31..7137df16b 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1367,7 +1367,7 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy if(alpha > 0.0){alpha /= compute_kinetic_energy(p_real, invM);} else{ alpha = 0.0; - cout << "Energy is drifted due to a classically forbidden hop" << endl; + cout << "Energy is drifted due to the dynamics initialization at a classical turning point" << endl; } for(int idof=0; idof Date: Wed, 20 Dec 2023 16:17:00 -0500 Subject: [PATCH 15/47] Added the function that returns the reprojection matrix based on state reordering and phase correction, as an alternative to local diabatization approach --- src/dyn/dyn_projectors.cpp | 49 +++++++++++++++++++++++++++++++++++++- src/dyn/dyn_projectors.h | 1 + 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/src/dyn/dyn_projectors.cpp b/src/dyn/dyn_projectors.cpp index 92afd9ffb..2988fce9b 100755 --- a/src/dyn/dyn_projectors.cpp +++ b/src/dyn/dyn_projectors.cpp @@ -1,5 +1,5 @@ /********************************************************************************* -* Copyright (C) 2019-2022 Alexey V. Akimov +* Copyright (C) 2019-2023 Alexey V. Akimov * * This file is distributed under the terms of the GNU General Public License * as published by the Free Software Foundation, either version 3 of @@ -761,6 +761,53 @@ vector compute_projectors(dyn_control_params& prms, vector& Ea } +CMATRIX compute_projector(dyn_control_params& prms, CMATRIX& Eadi, CMATRIX& St){ +/** + + Computes the instantaneous projector = permutation + phase correction + +*/ + + int nst = St.n_rows; + + vector perm_t(nst,0); + for(int i=0; i P * perm + res = permutation2cmatrix(perm_t); + st = st * res; + + if(prms.do_phase_correction){ + // ### Compute the instantaneous phase correction factors ### + phase_i = compute_phase_corrections(st, prms.phase_correction_tol); // f(i) + + // ### Scale projections' components by the phases ### + for(int a=0; a compute_projectors(dyn_control_params& prms, vector& St, vector >& perms){ /** diff --git a/src/dyn/dyn_projectors.h b/src/dyn/dyn_projectors.h index 8d91f79dc..0c3475f89 100755 --- a/src/dyn/dyn_projectors.h +++ b/src/dyn/dyn_projectors.h @@ -55,6 +55,7 @@ vector< vector > compute_permutations(dyn_control_params& prms, vector compute_projectors(dyn_control_params& prms, vector& Eadi, vector& St, Random& rnd); vector compute_projectors(dyn_control_params& prms, vector& St, vector >& perms); +CMATRIX compute_projector(dyn_control_params& prms, CMATRIX& Eadi, CMATRIX& St); CMATRIX raw_to_dynconsyst(CMATRIX& amplitudes, vector& projectors); CMATRIX dynconsyst_to_raw(CMATRIX& amplitudes, vector& projectors); From 7fbbdb714e450b1dd13318c3eeb8dbcb96c500e5 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Wed, 20 Dec 2023 16:17:58 -0500 Subject: [PATCH 16/47] exposed this function to Python --- src/dyn/libdyn.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index f80293338..bda60c87f 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -714,6 +714,9 @@ void export_dyn_projectors_objects(){ (dyn_control_params& prms, vector& St, vector >& perms) = &compute_projectors; def("compute_projectors", expt_compute_projectors_v2); + CMATRIX (*expt_compute_projector_v1) + (dyn_control_params& prms, CMATRIX& Eadi, CMATRIX& St) = &compute_projector; + def("compute_projector", expt_compute_projector_v1); CMATRIX (*expt_raw_to_dynconsyst_v1) (CMATRIX& amplitudes, vector& projectors) = &raw_to_dynconsyst; From c8de4365fdd4a63d0a91e07ff2544a542ee0ba8b Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Wed, 20 Dec 2023 16:19:39 -0500 Subject: [PATCH 17/47] Added the option to compute the basis re-projection matrix based on the state tracking + phase correction function in addition to the existing LD approach --- src/dyn/dyn_ham.cpp | 25 +++++++++++++++++++++---- 1 file changed, 21 insertions(+), 4 deletions(-) diff --git a/src/dyn/dyn_ham.cpp b/src/dyn/dyn_ham.cpp index 16cb972a3..bd535da0f 100644 --- a/src/dyn/dyn_ham.cpp +++ b/src/dyn/dyn_ham.cpp @@ -16,6 +16,7 @@ #include "dyn_ham.h" +#include "dyn_projectors.h" /// liblibra namespace namespace liblibra{ @@ -267,6 +268,8 @@ void update_proj_adi(dyn_control_params& prms, dyn_variables& dyn_var, nHamilton //======= Parameters of the dyn variables ========== int ntraj = dyn_var.ntraj; + double diff = 0.0; + for(int itraj=0; itraj=100 && method <200){ traj1 = 0; } @@ -279,14 +282,28 @@ void update_proj_adi(dyn_control_params& prms, dyn_variables& dyn_var, nHamilton P = ham->get_time_overlap_adi(); // (U_old.H() * U).H(); - - // More consistent with the new derivations: - FullPivLU_inverse(P, T_new); - T_new = orthogonalized_T( T_new ); + if(prms.state_tracking_algo==-1){ // This is LD + // More consistent with the new derivations: + FullPivLU_inverse(P, T_new); + + if( fabs( (P * T_new).tr().real() - P.n_cols) > 0.1 ){ + cout<<"Problem inverting time-overlap matrix\n"; + P.show_matrix("p_matrix.txt"); + exit(0); + } + T_new = orthogonalized_T( T_new ); + } + else{ // This is based on reordering + phase correction + CMATRIX Eadi(ham->get_ham_adi()); + T_new = P; + T_new = compute_projector(prms, Eadi, T_new); + } *dyn_var.proj_adi[itraj] = T_new; + }//for ntraj + }// reproject_basis From 1595e4ae30ee76ab3234d7d50282ce723349d528 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Wed, 20 Dec 2023 16:21:22 -0500 Subject: [PATCH 18/47] Redefined the state_tracking_algo to be -1 (local diabatization) as the default option --- src/dyn/dyn_control_params.cpp | 5 +++-- src/dyn/dyn_control_params.h | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/dyn/dyn_control_params.cpp b/src/dyn/dyn_control_params.cpp index 6078f3e3e..bd0dffc28 100755 --- a/src/dyn/dyn_control_params.cpp +++ b/src/dyn/dyn_control_params.cpp @@ -52,7 +52,7 @@ dyn_control_params::dyn_control_params(){ do_ssy = 0; do_phase_correction = 1; phase_correction_tol = 1e-3; - state_tracking_algo = 2; + state_tracking_algo = -1; MK_alpha = 0.0; MK_verbosity = 0; convergence = 0; @@ -217,7 +217,8 @@ dyn_control_params::~dyn_control_params() { void dyn_control_params::sanity_check(){ ///=================== Options for state tracking ====================== - if(state_tracking_algo==0 || state_tracking_algo==1 || + if(state_tracking_algo==-1 || + state_tracking_algo==0 || state_tracking_algo==1 || state_tracking_algo==2 || state_tracking_algo==3 || state_tracking_algo==32 || state_tracking_algo==33){ ; ; } else{ diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index b7ef90fe9..6d9163f28 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -244,9 +244,10 @@ class dyn_control_params{ /** State tracking algorithm: + - -1: use LD approach, it includes phase correction too [ default ] - 0: no state tracking - 1: method of Kosuke Sato (may fail by getting trapped into an infinite loop) - - 2: Munkres-Kuhn (Hungarian) algorithm [ default ] + - 2: Munkres-Kuhn (Hungarian) algorithm - 3: experimental stochastic algorithm, the original version with elimination (known problems) - 32: experimental stochastic algorithms with all permutations (too expensive) - 33: the improved stochastic algorithm with good scaling and performance, on par with the mincost From c7207368e17a2cc850f3416d3fbd997066c6c900 Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Wed, 20 Dec 2023 16:22:09 -0500 Subject: [PATCH 19/47] Added the default state_tracking_algo to be -1 to the Python code --- src/libra_py/dynamics/tsh/compute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index 750ddf9d3..583ecc61b 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -514,7 +514,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): "time_overlap_method":0, "nac_update_method":1, "nac_algo":0, "hvib_update_method":1, "do_ssy":0, "do_phase_correction":1, "phase_correction_tol":1e-3, - "state_tracking_algo":2, "MK_alpha":0.0, "MK_verbosity":0, + "state_tracking_algo":-1, "MK_alpha":0.0, "MK_verbosity":0, "convergence":0, "max_number_attempts":100, "min_probability_reordering":0.0, "is_nbra":0, "icond":0, "nfiles":-1, "thermally_corrected_nbra":0, "total_energy":0.01, "tcnbra_nu_therm":0.001, "tcnbra_nhc_size":1, "tcnbra_do_nac_scaling":1 From 5de4cf5a8193ee413f864c2f541f8ca26ad27739 Mon Sep 17 00:00:00 2001 From: MohammadShakiba Date: Sat, 23 Dec 2023 11:40:35 -0500 Subject: [PATCH 20/47] removed an auxiliary function from data_conv --- src/libra_py/data_conv.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/src/libra_py/data_conv.py b/src/libra_py/data_conv.py index 95dc85da1..34bc780cd 100755 --- a/src/libra_py/data_conv.py +++ b/src/libra_py/data_conv.py @@ -578,19 +578,4 @@ def vasp_to_xyz(filename): f.close() -def upper_vector_to_symmetric_nparray(upper_vector, upper_indices, mat_shape): - """ - This function gets the upper triangular part of a matrix as a vector and retuns a symmetric matrix - Args: - upper_vector (nparray): The upper triangular of a matrix - upper_indices (nparray): The indices of the upper part of the matrix - mat_shape (tuple): The shape of the original numpy array - Returns: - matrix (nparray): The symmetric matix built based on the upper triangular matrix - """ - matrix = np.zeros(mat_shape) - matrix[upper_indices] = upper_vector - matrix = matrix + matrix.T - np.diag(matrix.diagonal()) - return matrix - From 639da5b6651ab67701b110b4fdaac5af02360125 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sat, 30 Dec 2023 03:14:17 -0500 Subject: [PATCH 21/47] Add more turning-point algorithms --- src/dyn/dyn_control_params.cpp | 3 + src/dyn/dyn_control_params.h | 10 ++ src/dyn/dyn_decoherence_methods.cpp | 149 +++++++++++++++++++++------ src/dyn/dyn_variables.cpp | 8 ++ src/dyn/dyn_variables.h | 16 +++ src/dyn/libdyn.cpp | 1 + src/libra_py/dynamics/tsh/compute.py | 8 +- 7 files changed, 162 insertions(+), 33 deletions(-) diff --git a/src/dyn/dyn_control_params.cpp b/src/dyn/dyn_control_params.cpp index bd0dffc28..00ea916a4 100755 --- a/src/dyn/dyn_control_params.cpp +++ b/src/dyn/dyn_control_params.cpp @@ -83,6 +83,7 @@ dyn_control_params::dyn_control_params(){ coherence_threshold = 0.01; use_xf_force = 0; project_out_aux = 0; + tp_algo = 1; ///================= Entanglement of trajectories ================================ entanglement_opt = 0; @@ -163,6 +164,7 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){ coherence_threshold = x.coherence_threshold; use_xf_force = x.use_xf_force; project_out_aux = x.project_out_aux; + tp_algo = x.tp_algo; ///================= Entanglement of trajectories ================================ entanglement_opt = x.entanglement_opt; @@ -350,6 +352,7 @@ void dyn_control_params::set_parameters(bp::dict params){ else if(key=="coherence_threshold"){ coherence_threshold = bp::extract(params.values()[i]); } else if(key=="use_xf_force"){ use_xf_force = bp::extract(params.values()[i]); } else if(key=="project_out_aux"){ project_out_aux = bp::extract(params.values()[i]); } + else if(key=="tp_algo"){ tp_algo = bp::extract(params.values()[i]); } ///================= Entanglement of trajectories ================================ else if(key=="entanglement_opt"){ entanglement_opt = bp::extract(params.values()[i]); } diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index 6d9163f28..9a19c716b 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -528,6 +528,16 @@ class dyn_control_params{ */ int project_out_aux; + /** + Turning-point algorithm for auxiliary trajectories + + Options: + - 0: no treatment of a turning point + - 1: fix auxiliary trajectories of adiabatic states other than the active state [default] + - 2: collapse to the active state + */ + int tp_algo; + /** 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. diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 7137df16b..2289ee679 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -949,12 +949,16 @@ void xf_init_AT(dyn_variables& dyn_var, int traj, int ist){ if( ist == -1){ dyn_var.is_mixed[traj].assign(nadi, 0); dyn_var.is_first[traj].assign(nadi, 0); + dyn_var.is_fixed[traj].assign(nadi, 0); + dyn_var.is_keep[traj].assign(nadi, 0); dyn_var.p_quant->set(-1, traj, 0.0); dyn_var.VP->set(-1, traj, 0.0); } else{ dyn_var.is_mixed[traj][ist] = 0; dyn_var.is_first[traj][ist] = 0; + dyn_var.is_fixed[traj][ist] = 0; + dyn_var.is_keep[traj][ist] = 0; } } @@ -1118,6 +1122,7 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn vector& is_mixed = dyn_var.is_mixed[traj]; vector& is_first = dyn_var.is_first[traj]; + vector& is_fixed = dyn_var.is_fixed[traj]; MATRIX& q_aux = *dyn_var.q_aux[traj]; MATRIX& p_aux = *dyn_var.p_aux[traj]; @@ -1126,6 +1131,8 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn for(int i=0; iget(idof, traj));} @@ -1148,6 +1155,8 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; vector& is_first = dyn_var.is_first[traj]; + vector& is_fixed = dyn_var.is_fixed[traj]; + vector& is_keep = dyn_var.is_keep[traj]; MATRIX& p_aux = *dyn_var.p_aux[traj]; MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; @@ -1167,9 +1176,12 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn } } + int is_tp = 0; double alpha; for(int i=0; icol(traj); if(i==a){ @@ -1200,19 +1212,50 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn } // Check the turning point - if (is_first[i] == 0){ + if (is_first[i] == 0 and prms.tp_algo != 0){ double temp = 0.0; for(int idof=0; idof& is_mixed = dyn_var.is_mixed[traj]; vector& is_first = dyn_var.is_first[traj]; + vector& is_fixed = dyn_var.is_fixed[traj]; MATRIX& q_aux = *dyn_var.q_aux[traj]; MATRIX& p_aux = *dyn_var.p_aux[traj]; for(int i=0; iget(idof, traj));} @@ -1285,6 +1331,8 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; vector& is_first = dyn_var.is_first[traj]; + vector& is_fixed = dyn_var.is_fixed[traj]; + vector& is_keep = dyn_var.is_keep[traj]; MATRIX& p_aux = *dyn_var.p_aux[traj]; MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; @@ -1305,10 +1353,13 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy vector _id(2, 0); _id[1] = traj; coeff = dyn_var.ampl_adi->col(traj); double Epot = ham.Ehrenfest_energy_adi(coeff, _id).real(); - + + int is_tp = 0; double alpha; for(int i=0; icol(traj); if(is_first[i]==1){ @@ -1347,41 +1398,75 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy } // Check the turning point - if (is_first[i] == 0){ + if (is_first[i] == 0 and prms.tp_algo != 0){ double temp = 0.0; for(int idof=0; idofcol(traj); - double Epot = ham.Ehrenfest_energy_adi(coeff, _id).real(); - - // Rescaling momenta for the energy conservation - p_real = dyn_var.p->col(traj); - double alpha = compute_kinetic_energy(p_real, invM) + Epot_old - Epot; + } + } - if(alpha > 0.0){alpha /= compute_kinetic_energy(p_real, invM);} - else{ - alpha = 0.0; - cout << "Energy is drifted due to the dynamics initialization at a classical turning point" << endl; - } + } + }//i + + if(is_tp == 1){ + if(prms.tp_algo == 1){ + int a = dyn_var.act_states[traj]; + for(int i=0; iset(idof, traj, dyn_var.p->get(idof, traj) * sqrt(alpha)); - } + int a = dyn_var.act_states[traj]; + collapse(*dyn_var.ampl_adi, traj, a, 0); + + // After the collapse + coeff = dyn_var.ampl_adi->col(traj); + double Epot = ham.Ehrenfest_energy_adi(coeff, _id).real(); + + // Rescaling momenta for the energy conservation + p_real = dyn_var.p->col(traj); + double alpha = compute_kinetic_energy(p_real, invM) + Epot_old - Epot; - xf_init_AT(dyn_var, traj, -1); - cout << "Collapse to the active state " << a << " at a classical turning point on traj " << traj < 0.0){alpha /= compute_kinetic_energy(p_real, invM);} + else{ + alpha = 0.0; + cout << "Energy is drifted due to the dynamics initialization at a classical turning point" << endl; + } + + for(int idof=0; idofset(idof, traj, dyn_var.p->get(idof, traj) * sqrt(alpha)); } + xf_init_AT(dyn_var, traj, -1); + cout << "Collapse to the active state " << a << " at a classical turning point on traj " << traj <col(traj); double Ekin = compute_kinetic_energy(p_real, invM); - double e = 1.0e-6; // masking for classical turning points + double e = 1.0e-4; // masking for classical turning points for(int i=0; i()); is_first.push_back(vector()); + is_fixed.push_back(vector()); + is_keep.push_back(vector()); for(int i=0; i()); is_first.push_back(vector()); + is_fixed.push_back(vector()); + is_keep.push_back(vector()); for(int i=0; i > is_first(ntraj, nadi) */ vector> is_first; + + /** + Whether to fix an auxiliary trajectory + + Options: + vector< vector > is_fixed(ntraj, nadi) + */ + vector> is_fixed; + + /** + Whether to keep the auxiliary momenta + + Options: + vector< vector > is_keep(ntraj, nadi) + */ + vector> is_keep; /** Nuclear coordinates of state-wise auxiliary trajectories diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index bda60c87f..0cd63cf29 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -105,6 +105,7 @@ void export_dyn_control_params_objects(){ .def_readwrite("coherence_threshold", &dyn_control_params::coherence_threshold) .def_readwrite("use_xf_force", &dyn_control_params::use_xf_force) .def_readwrite("project_out_aux", &dyn_control_params::project_out_aux) + .def_readwrite("tp_algo", &dyn_control_params::tp_algo) ///================= Entanglement of trajectories ================================ .def_readwrite("entanglement_opt", &dyn_control_params::entanglement_opt) diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index 583ecc61b..d6e54d628 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -322,6 +322,12 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): * **dyn_params["project_out_aux"]** (int): Whether to project out the density on an auxiliary trajectory when its motion is classically forbidden. [ default: 0] Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` + * **dyn_params["tp_algo"]** (int): Turning-point algorithm for auxiliary trajectories + + - 0: no treatment of a turning point + - 1: fix auxiliary trajectories of adiabatic states other than the active state [default] + - 2: collapse to the active state + ///=============================================================================== ///================= Entanglement of trajectories ================================ ///=============================================================================== @@ -535,7 +541,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): "ave_gaps":MATRIX(nstates,nstates), "schwartz_decoherence_inv_alpha": MATRIX(nstates, 1), "wp_width":0.3, "coherence_threshold":0.01, "use_xf_force": 0, - "project_out_aux": 0 + "project_out_aux": 0, "tp_algo": 1 } ) #================= Entanglement of trajectories ================================ From ac23286474487978d4d797d2ab96fb4e41938a28 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sun, 31 Dec 2023 16:56:25 -0500 Subject: [PATCH 22/47] Organize turning-point algorithms --- src/dyn/dyn_control_params.h | 5 ++-- src/dyn/dyn_decoherence_methods.cpp | 36 ++++++++++++++-------------- src/libra_py/dynamics/tsh/compute.py | 5 ++-- 3 files changed, 24 insertions(+), 22 deletions(-) diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index 9a19c716b..ad1b6ff37 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -533,8 +533,9 @@ class dyn_control_params{ Options: - 0: no treatment of a turning point - - 1: fix auxiliary trajectories of adiabatic states other than the active state [default] - - 2: collapse to the active state + - 1: collapse to the active state [default] + - 2: fix auxiliary positions of adiabatic states except for the active state + - 3: keep auxiliary momenta of adiabatic states except for the active state */ int tp_algo; diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 2289ee679..becbb255c 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1226,6 +1226,11 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn if(is_tp == 1){ if(prms.tp_algo == 1){ + collapse(*dyn_var.ampl_adi, traj, a, 0); + xf_init_AT(dyn_var, traj, -1); + cout << "Collapse to the active state " << a << " at a classical turning point on traj " << traj < Date: Sun, 31 Dec 2023 19:19:33 -0500 Subject: [PATCH 23/47] Simplify the XF propagation --- src/dyn/Dynamics.cpp | 4 ++-- src/dyn/dyn_decoherence.h | 2 +- src/dyn/dyn_decoherence_methods.cpp | 22 ++++++++++++++-------- 3 files changed, 17 insertions(+), 11 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 6fc13cd8a..5cef3211b 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1217,13 +1217,13 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, // basis re-projection matrices for(i=0; i 0.0){alpha /= compute_kinetic_energy(p_real, *dyn_var.iM);} - else{alpha = 0.0;} + else{alpha = 0.0; + cout << "Total energy is drifted due to adiabatic recovery to a classically forbidden state!" << endl; + } for(int idof=0; idofset(idof, traj, dyn_var.p->get(idof, traj) * sqrt(alpha)); @@ -1429,7 +1431,7 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy if(alpha > 0.0){alpha /= compute_kinetic_energy(p_real, invM);} else{ alpha = 0.0; - cout << "Energy is drifted due to the dynamics initialization at a classical turning point" << endl; + cout << "Total energy is drifted due to the dynamics initialization at a classical turning point" << endl; } for(int idof=0; idof(0.0, -0.25*dt) ); + D = libspecialfunctions::exp_(Hxf, complex(0.0, -0.5*dt) ); + if(rotation == 1){ + D = T_new * D * T_new.H(); + } C = D * C; From 82b69fa98f38bcf7e2573a2e083d2fe4f35bff35 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sun, 31 Dec 2023 20:03:37 -0500 Subject: [PATCH 24/47] Add the masking constant to param --- src/dyn/dyn_control_params.cpp | 3 +++ src/dyn/dyn_control_params.h | 13 ++++++++++++- src/dyn/dyn_decoherence_methods.cpp | 3 +-- src/dyn/libdyn.cpp | 1 + src/libra_py/dynamics/tsh/compute.py | 8 ++++++-- 5 files changed, 23 insertions(+), 5 deletions(-) diff --git a/src/dyn/dyn_control_params.cpp b/src/dyn/dyn_control_params.cpp index 00ea916a4..aa9cd21f1 100755 --- a/src/dyn/dyn_control_params.cpp +++ b/src/dyn/dyn_control_params.cpp @@ -81,6 +81,7 @@ dyn_control_params::dyn_control_params(){ ave_gaps = NULL; wp_width = 0.1; coherence_threshold = 0.01; + e_mask = 0.0001; use_xf_force = 0; project_out_aux = 0; tp_algo = 1; @@ -162,6 +163,7 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){ collapse_option = x.collapse_option; wp_width = x.wp_width; coherence_threshold = x.coherence_threshold; + e_mask = x.e_mask; use_xf_force = x.use_xf_force; project_out_aux = x.project_out_aux; tp_algo = x.tp_algo; @@ -350,6 +352,7 @@ void dyn_control_params::set_parameters(bp::dict params){ } else if(key=="wp_width"){ wp_width = bp::extract(params.values()[i]); } else if(key=="coherence_threshold"){ coherence_threshold = bp::extract(params.values()[i]); } + else if(key=="e_mask"){ e_mask = bp::extract(params.values()[i]); } else if(key=="use_xf_force"){ use_xf_force = bp::extract(params.values()[i]); } else if(key=="project_out_aux"){ project_out_aux = bp::extract(params.values()[i]); } else if(key=="tp_algo"){ tp_algo = bp::extract(params.values()[i]); } diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index ad1b6ff37..e49f2aa3e 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -514,7 +514,15 @@ class dyn_control_params{ [ default : 0.01 ] */ double coherence_threshold; - + + + /** + The masking parameter for computing nabla phase vectors in the MQCXF + [ default : 0.0001 ] + */ + double e_mask; + + /** Whether to use the decoherence force in MQCXF The corresponding electronic propagation is adjusted for the energy conservation @@ -522,12 +530,14 @@ class dyn_control_params{ */ int use_xf_force; + /** Whether to project out the density on an auxiliary trajectory when its motion is classically forbidden [ default : 0 ] */ int project_out_aux; + /** Turning-point algorithm for auxiliary trajectories @@ -539,6 +549,7 @@ class dyn_control_params{ */ int tp_algo; + /** 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. diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index ccd13467f..12ed6c77b 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1484,11 +1484,10 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy p_real = dyn_var.p->col(traj); double Ekin = compute_kinetic_energy(p_real, invM); - double e = 1.0e-4; // masking for classical turning points for(int i=0; iget(idof,traj)/(Ekin + e)); + nab_phase.set(i, idof, -0.5*E.get(i,i).real()*dyn_var.p->get(idof,traj)/(Ekin + prms.e_mask)); }//idof } }//i diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index 0cd63cf29..941deac59 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -103,6 +103,7 @@ void export_dyn_control_params_objects(){ .def_readwrite("ave_gaps", &dyn_control_params::ave_gaps) .def_readwrite("wp_width", &dyn_control_params::wp_width) .def_readwrite("coherence_threshold", &dyn_control_params::coherence_threshold) + .def_readwrite("e_mask", &dyn_control_params::e_mask) .def_readwrite("use_xf_force", &dyn_control_params::use_xf_force) .def_readwrite("project_out_aux", &dyn_control_params::project_out_aux) .def_readwrite("tp_algo", &dyn_control_params::tp_algo) diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index 4c077bd4a..c3c81156b 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -313,6 +313,10 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` + * **dyn_params["e_mask"]** ( double ): The masking parameter for computing nabla phase vectors. [ default: 0.0001 ] + Only used with the MQCXF method, that is, `decoherence_algo == 5` + + * **dyn_params["use_xf_force"]** (int): Whether to use the XF-based force. Only used with `decoherence_algo == 6` @@ -541,8 +545,8 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): "decoherence_rates":MATRIX(nstates, nstates), "ave_gaps":MATRIX(nstates,nstates), "schwartz_decoherence_inv_alpha": MATRIX(nstates, 1), - "wp_width":0.3, "coherence_threshold":0.01, "use_xf_force": 0, - "project_out_aux": 0, "tp_algo": 1 + "wp_width":0.3, "coherence_threshold":0.01, "e_mask": 0.0001, + "use_xf_force": 0, "project_out_aux": 0, "tp_algo": 1 } ) #================= Entanglement of trajectories ================================ From 45c3bb207a304cfbe055236aa5c02a22197d00b6 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sun, 31 Dec 2023 20:54:02 -0500 Subject: [PATCH 25/47] Modify p_aux computation due to total energy jumps --- src/dyn/dyn_decoherence_methods.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 12ed6c77b..35da6e904 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1368,9 +1368,9 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy alpha = compute_kinetic_energy(p_real, invM) + Epot - ham_adi.get(i,i).real(); } else{ - p_aux_temp = p_aux_old.row(i).T(); - alpha = compute_kinetic_energy(p_aux_temp, invM) + ham_adi_prev.get(i,i).real() - ham_adi.get(i,i).real(); - //alpha = compute_kinetic_energy(p_real, invM) + Epot - ham_adi.get(i,i).real(); + //p_aux_temp = p_aux_old.row(i).T(); + //alpha = compute_kinetic_energy(p_aux_temp, invM) + ham_adi_prev.get(i,i).real() - ham_adi.get(i,i).real(); + alpha = compute_kinetic_energy(p_real, invM) + Epot - ham_adi.get(i,i).real(); } if (alpha < 0.0){ alpha = 0.0; From bbb40900d8872fb3560d7c45c30993ebefa0c4e5 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sun, 31 Dec 2023 23:12:31 -0500 Subject: [PATCH 26/47] Add missing initialization of auxiliary trajectories --- src/dyn/dyn_decoherence_methods.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 35da6e904..01024c58c 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1075,7 +1075,10 @@ void xf_create_AT(dyn_variables& dyn_var, double threshold){ int nr_mixed = 0; for(int i=0; iupper_lim){is_mixed[i]=0;} + if(a_ii<=lower_lim || a_ii>upper_lim){ + is_mixed[i]=0; + xf_init_AT(dyn_var, traj, i); + } else{ is_mixed[i]==1 ? is_first[i]=0:is_first[i]=1; is_mixed[i]=1; From ab4c54909d98e06a91ac52a190f4e0e8e825d7c1 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Mon, 1 Jan 2024 01:17:04 -0500 Subject: [PATCH 27/47] Add safety related to is_mixed --- src/dyn/dyn_decoherence_methods.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 01024c58c..1b54212de 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1529,7 +1529,9 @@ void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_w // Set a diagonal matrix of nabla_phase for each dof CMATRIX F(nst, nst); for(int i=0; i(0.0, dyn_var.nab_phase[traj]->get(i, idof)) ); + if(is_mixed[i]==1){ + F.set(i,i, complex(0.0, dyn_var.nab_phase[traj]->get(i, idof)) ); + } } F = T * F * T.H(); @@ -1573,10 +1575,13 @@ void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& h double Ekin = compute_kinetic_energy(p_real, invM); // Compute F for each dof + vector& is_mixed = dyn_var.is_mixed[traj]; for(int idof=0; idof (0.0, 0.0)); for(int i=0; i (dyn_var.nab_phase[traj]->get(i, idof), 0.0) ); + if(is_mixed[i]==1){ + F[idof].set(i,i,complex (dyn_var.nab_phase[traj]->get(i, idof), 0.0) ); + } } } From 31ee07da9cfbddbe711517974df543a38ba306aa Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Mon, 1 Jan 2024 06:08:51 -0500 Subject: [PATCH 28/47] Add a crude td wavepacket width functionality --- src/dyn/Dynamics.cpp | 4 ++-- src/dyn/dyn_control_params.cpp | 6 ++++++ src/dyn/dyn_control_params.h | 20 ++++++++++++++++++- src/dyn/dyn_decoherence.h | 2 +- src/dyn/dyn_decoherence_methods.cpp | 30 ++++++++++++++++++++++++++-- src/libra_py/dynamics/tsh/compute.py | 22 +++++++++++++++++--- 6 files changed, 75 insertions(+), 9 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 5cef3211b..c5838f4cc 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1217,13 +1217,13 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, // basis re-projection matrices for(i=0; i(params.values()[i]); } + else if(key=="wp_v"){ wp_v = bp::extract(params.values()[i]); } else if(key=="coherence_threshold"){ coherence_threshold = bp::extract(params.values()[i]); } else if(key=="e_mask"){ e_mask = bp::extract(params.values()[i]); } else if(key=="use_xf_force"){ use_xf_force = bp::extract(params.values()[i]); } else if(key=="project_out_aux"){ project_out_aux = bp::extract(params.values()[i]); } else if(key=="tp_algo"){ tp_algo = bp::extract(params.values()[i]); } + else if(key=="use_td_width"){ use_td_width = bp::extract(params.values()[i]); } ///================= Entanglement of trajectories ================================ else if(key=="entanglement_opt"){ entanglement_opt = bp::extract(params.values()[i]); } diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index e49f2aa3e..74719a17b 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -503,12 +503,20 @@ class dyn_control_params{ /** - The width of frozen Gaussian for the decoherence from SHXF & MQCXF + The width of Gaussian for the decoherence from SHXF & MQCXF [ default : 0.1 ] */ double wp_width; + /** + The velocity of Gaussian for the decoherence from SHXF & MQCXF + This value is applied when use_td_width = 1 + [ default : 0.0 ] + */ + double wp_v; + + /** The criterion whether the electronic state is in a coherence [ default : 0.01 ] @@ -548,6 +556,16 @@ class dyn_control_params{ - 3: keep auxiliary momenta of adiabatic states except for the active state */ int tp_algo; + + + /** + Whether to use the td Gaussian width for the nuclear wave packet approximation + This option can be considered when it comes to unbounded systems. + This approximation is based on a nuclear wave packet on a free surface: + \sigma_x(t)=\sqrt[\sigma_x(0)^2 + (wp_v * t)^2] + [ default : 0 ] + */ + int use_td_width; /** diff --git a/src/dyn/dyn_decoherence.h b/src/dyn/dyn_decoherence.h index 35b9bebe3..dee81c288 100755 --- a/src/dyn/dyn_decoherence.h +++ b/src/dyn/dyn_decoherence.h @@ -81,7 +81,7 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy // XF propagation void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev); -void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms, int rotation); +void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms, bp::dict params, int rotation); void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj); ///================ In dyn_decoherence_time.cpp =================================== diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 1b54212de..4ada59eef 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1497,6 +1497,24 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy } // traj } +double compute_wp_width(dyn_variables& dyn_var, dyn_control_params& prms, bp::dict params, double dt){ + double elapsed_time, s2; + int it; + + std::string key; + for(int i=0;i(params.keys()[i]); + if(key=="timestep") { it = bp::extract(params.values()[i]); } + else {continue;} + } + + elapsed_time = dt*it; + s2 = pow(prms.wp_width, 2.0) + pow(prms.wp_v*elapsed_time, 2.0); + + return sqrt(s2); + +} + void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj){ int ndof = dyn_var.ndof; @@ -1615,7 +1633,7 @@ void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& h *dyn_var.f += *dyn_var.f_xf; } -void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms, int rotation){ +void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms, bp::dict params, int rotation){ int itraj, i, j; int num_el = prms.num_electronic_substeps; @@ -1653,7 +1671,15 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa CMATRIX Hxf(nadi, nadi); CMATRIX D(nadi, nadi); /// this is \exp[-idt/4\hbar * ( T_new.H()*Hxf(t+dt)*T_new + Hxf(t) )] - XF_correction(Hxf_old, dyn_var, C, prms.wp_width, T, itraj); + double wp_width; + if (prms.use_td_width == 1){ + wp_width = compute_wp_width(dyn_var, prms, params, dt); + } + else{ + wp_width = prms.wp_width; + } + + XF_correction(Hxf_old, dyn_var, C, wp_width, T, itraj); //XF_correction(Hxf, dyn_var, C, prms.wp_width, T, itraj); //Hxf = T_new.H() * Hxf * T_new; diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index c3c81156b..5d02f42bd 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -306,6 +306,11 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): * **dyn_params["wp_width"]** ( double ): A width of a Gaussian function as an approximation to adiabatic wave packets. [ default: 0.3 Bohr ] + This value is used for the initial conditions when the td Gaussian approximation is used, that is, use_td_width == 1 + Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` + + + * **dyn_params["wp_v"]** ( double ): The velocity of Gaussian wave packet following potential-free td Gaussian [ default: 0.0 a.u. ] Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` @@ -313,7 +318,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` - * **dyn_params["e_mask"]** ( double ): The masking parameter for computing nabla phase vectors. [ default: 0.0001 ] + * **dyn_params["e_mask"]** ( double ): The masking parameter for computing nabla phase vectors. [ default: 0.0001 Ha ] Only used with the MQCXF method, that is, `decoherence_algo == 5` @@ -323,16 +328,27 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): - 0: Only Ehrenfest-like force; EhXF [ default ] - 1: The whole force including the XF-correction; MQCXF + * **dyn_params["project_out_aux"]** (int): Whether to project out the density on an auxiliary trajectory when its motion is classically forbidden. [ default: 0] Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` + * **dyn_params["tp_algo"]** (int): Turning-point algorithm for auxiliary trajectories + Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` - 0: no treatment of a turning point - 1: collapse to the active state [default] - 2: fix auxiliary positions of adiabatic states except for the active state - 3: keep auxiliary momenta of adiabatic states except for the active state + + * **dyn_params["use_td_width"]** (int): Whether to use the td Gaussian width for the nuclear wave packet approximation [ default : 0 ] + This option can be considered when it comes to unbounded systems. + This approximation is based on a nuclear wave packet on a free surface: + \sigma_x(t)=\sqrt[\sigma_x(0)^2 + (wp_v * t)^2] + Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` + + ///=============================================================================== ///================= Entanglement of trajectories ================================ ///=============================================================================== @@ -545,8 +561,8 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): "decoherence_rates":MATRIX(nstates, nstates), "ave_gaps":MATRIX(nstates,nstates), "schwartz_decoherence_inv_alpha": MATRIX(nstates, 1), - "wp_width":0.3, "coherence_threshold":0.01, "e_mask": 0.0001, - "use_xf_force": 0, "project_out_aux": 0, "tp_algo": 1 + "wp_width":0.3, "wp_v":0, "coherence_threshold":0.01, "e_mask": 0.0001, + "use_xf_force": 0, "project_out_aux": 0, "tp_algo": 1, "use_td_width": 0 } ) #================= Entanglement of trajectories ================================ From 3115877907dab2ce1cca6e49cad314e34c20ff31 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Wed, 3 Jan 2024 00:04:52 -0500 Subject: [PATCH 29/47] Add a function to compute the coherence in the exact calculation --- src/dyn/wfcgrid2/Wfcgrid2.h | 2 + src/dyn/wfcgrid2/Wfcgrid2_properties.cpp | 59 ++++++++++++++++++++++++ src/dyn/wfcgrid2/libwfcgrid2.cpp | 4 ++ src/libra_py/dynamics/exact/save.py | 8 ++++ 4 files changed, 73 insertions(+) diff --git a/src/dyn/wfcgrid2/Wfcgrid2.h b/src/dyn/wfcgrid2/Wfcgrid2.h index 706c6a45f..140cf9b80 100755 --- a/src/dyn/wfcgrid2/Wfcgrid2.h +++ b/src/dyn/wfcgrid2/Wfcgrid2.h @@ -223,6 +223,8 @@ class Wfcgrid2{ MATRIX get_pops(int rep); MATRIX get_pops(int rep, vector& bmin, vector& bmax); + + MATRIX get_coherences(int rep); void compute_wfc_gradients(int rep, int idof, double mass); diff --git a/src/dyn/wfcgrid2/Wfcgrid2_properties.cpp b/src/dyn/wfcgrid2/Wfcgrid2_properties.cpp index 4191dff4c..e55e67ae2 100755 --- a/src/dyn/wfcgrid2/Wfcgrid2_properties.cpp +++ b/src/dyn/wfcgrid2/Wfcgrid2_properties.cpp @@ -423,6 +423,65 @@ MATRIX Wfcgrid2::get_pops(int rep, vector& bmin, vector& bmax){ } +MATRIX Wfcgrid2::get_coherences(int rep){ +/** + Compute the coherence indicator between all states in a given representation + + Out: MATRIX(nstates, nstates) + +*/ + MATRIX res(nstates, nstates); + + int istate, jstate, npt1; + + MATRIX temp(Npts, 1); + if(rep=0){ + for(npt1=0; npt1& bmin, vector& bmax) = &Wfcgrid2::get_pops; + MATRIX (Wfcgrid2::*expt_get_coherences_v1) + (int rep) = &Wfcgrid2::get_coherences; + // Auxiliary functions @@ -217,6 +220,7 @@ void export_Wfcgrid2_objects(){ .def("get_den_mat", &Wfcgrid2::get_den_mat) .def("get_pops", expt_get_pops_v1) .def("get_pops", expt_get_pops_v2) + .def("get_coherences", expt_get_coherences_v1) /** Wfcgrid2_SOFT */ .def("update_propagator_H", &Wfcgrid2::update_propagator_H) diff --git a/src/libra_py/dynamics/exact/save.py b/src/libra_py/dynamics/exact/save.py index 9923df513..0eb43522c 100755 --- a/src/libra_py/dynamics/exact/save.py +++ b/src/libra_py/dynamics/exact/save.py @@ -129,6 +129,12 @@ def exact_init_hdf5(saver, hdf5_output_level, _nsteps, _ndof, _nstates, _ngrid): # Density matrix computed in adiabatic rep saver.add_dataset("denmat_adi", (_nsteps, _nstates, _nstates), "C") + + # Density matrix computed in adiabatic rep + saver.add_dataset("coherence_dia", (_nsteps, _nstates, _nstates), "R") + + # Density matrix computed in adiabatic rep + saver.add_dataset("coherence_adi", (_nsteps, _nstates, _nstates), "R") if hdf5_output_level>=4: @@ -201,6 +207,8 @@ def save_data_hdf5(step, wfc, saver, params): if hdf5_output_level>=3: saver.save_matrix(step, "denmat_dia", wfc.get_den_mat(0) ) saver.save_matrix(step, "denmat_adi", wfc.get_den_mat(1) ) + saver.save_matrix(step, "coherence_dia", wfc.get_coherences(0) ) + saver.save_matrix(step, "coherence_adi", wfc.get_coherences(1) ) if hdf5_output_level>=4: From 2c7e0224f6281aaad8b0d4987fa627296d9b8ac5 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Wed, 3 Jan 2024 20:31:59 -0500 Subject: [PATCH 30/47] Make separate variables for current timestep and wp_width and related functions --- src/dyn/Dynamics.cpp | 39 ++++++++++++++++++++++++++-- src/dyn/dyn_control_params.cpp | 29 ++++++++++++++++----- src/dyn/dyn_control_params.h | 12 ++++----- src/dyn/dyn_decoherence.h | 2 +- src/dyn/dyn_decoherence_methods.cpp | 38 ++++----------------------- src/dyn/dyn_variables.cpp | 6 +++++ src/dyn/dyn_variables.h | 25 +++++++++++++++++- src/dyn/libdyn.cpp | 1 + src/libra_py/dynamics/tsh/compute.py | 6 ++--- 9 files changed, 106 insertions(+), 52 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index c5838f4cc..5a867eac0 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -570,6 +570,36 @@ void remove_thermal_correction(dyn_variables& dyn_var, nHamiltonian& ham, dyn_co } +void update_wp_width(dyn_variables& dyn_var, dyn_control_params& prms){ +/** + Updates the wave packet width in XF based on the Gaussian approximation +*/ + + //======= Parameters of the dyn variables ========== + int ndof = dyn_var.ndof; + int ntraj = dyn_var.ntraj; + + if (prms.use_td_width == 0){ + for(int itraj=0; itrajset(idof, itraj, prms.wp_width->get(idof, 0)); + } + } + } + else if (prms.use_td_width == 1){ + double elapsed_time, s2; + + elapsed_time = prms.dt*dyn_var.timestep; + + for(int itraj=0; itrajget(idof, 0), 2.0) + pow(prms.wp_v->get(idof, 0)*elapsed_time, 2.0); + dyn_var.wp_width->set(idof, itraj, sqrt(s2)); + } + } + } +} + /* void update_proj_adi(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonian* Ham_prev, dyn_control_params& prms){ @@ -1212,18 +1242,23 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, // Recompute NAC, Hvib, etc. in response to change of p update_Hamiltonian_variables(prms, dyn_var, ham, ham_aux, py_funct, params, 1); + // Update wave packet width in XF algorithm + if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ + dyn_var.get_current_timestep(params); + update_wp_width(dyn_var, prms); + } // Propagate electronic coefficients in the [t, t + dt] interval, this also updates the // basis re-projection matrices for(i=0; in_rows, x.wp_width->n_cols ); + *wp_width = *x.wp_width; + wp_v = new MATRIX( x.wp_v->n_rows, x.wp_v->n_cols ); + *wp_v = *x.wp_v; coherence_threshold = x.coherence_threshold; e_mask = x.e_mask; use_xf_force = x.use_xf_force; @@ -216,6 +219,8 @@ dyn_control_params::~dyn_control_params() { //cout<<"dyn_control_params destructor\n"; + delete wp_width; + delete wp_v; delete decoherence_rates; delete ave_gaps; delete schwartz_decoherence_inv_alpha; @@ -354,8 +359,20 @@ void dyn_control_params::set_parameters(bp::dict params){ for(int b=0;bset(a, b, x.get(a,b)); } } } - else if(key=="wp_width"){ wp_width = bp::extract(params.values()[i]); } - else if(key=="wp_v"){ wp_v = bp::extract(params.values()[i]); } + else if(key=="wp_width"){ + MATRIX x( bp::extract(params.values()[i]) ); + wp_width = new MATRIX(x.n_rows, x.n_cols); + for(int a=0;aset(a, b, x.get(a,b)); } + } + } + else if(key=="wp_v"){ + MATRIX x( bp::extract(params.values()[i]) ); + wp_v = new MATRIX(x.n_rows, x.n_cols); + for(int a=0;aset(a, b, x.get(a,b)); } + } + } else if(key=="coherence_threshold"){ coherence_threshold = bp::extract(params.values()[i]); } else if(key=="e_mask"){ e_mask = bp::extract(params.values()[i]); } else if(key=="use_xf_force"){ use_xf_force = bp::extract(params.values()[i]); } diff --git a/src/dyn/dyn_control_params.h b/src/dyn/dyn_control_params.h index 74719a17b..90152b9a1 100755 --- a/src/dyn/dyn_control_params.h +++ b/src/dyn/dyn_control_params.h @@ -503,18 +503,18 @@ class dyn_control_params{ /** - The width of Gaussian for the decoherence from SHXF & MQCXF - [ default : 0.1 ] + MATRIX(ndof, 1) of (initial) wave packet widths for the decoherence from SHXF & MQCXF + [ default : NULL ] */ - double wp_width; + MATRIX* wp_width; /** - The velocity of Gaussian for the decoherence from SHXF & MQCXF + MATRIX(ndof, 1) of wave packet velocities for the decoherence from SHXF & MQCXF This value is applied when use_td_width = 1 - [ default : 0.0 ] + [ default : NULL ] */ - double wp_v; + MATRIX* wp_v; /** diff --git a/src/dyn/dyn_decoherence.h b/src/dyn/dyn_decoherence.h index dee81c288..35b9bebe3 100755 --- a/src/dyn/dyn_decoherence.h +++ b/src/dyn/dyn_decoherence.h @@ -81,7 +81,7 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy // XF propagation void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev); -void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms, bp::dict params, int rotation); +void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms, int rotation); void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj); ///================ In dyn_decoherence_time.cpp =================================== diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 4ada59eef..dd77c49e9 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1497,25 +1497,7 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy } // traj } -double compute_wp_width(dyn_variables& dyn_var, dyn_control_params& prms, bp::dict params, double dt){ - double elapsed_time, s2; - int it; - - std::string key; - for(int i=0;i(params.keys()[i]); - if(key=="timestep") { it = bp::extract(params.values()[i]); } - else {continue;} - } - - elapsed_time = dt*it; - s2 = pow(prms.wp_width, 2.0) + pow(prms.wp_v*elapsed_time, 2.0); - - return sqrt(s2); - -} - -void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj){ +void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, CMATRIX& T, int traj){ int ndof = dyn_var.ndof; int nst = dyn_var.nadi; @@ -1534,7 +1516,7 @@ void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_w for(int i=0; iadd(idof, traj, 0.5 / pow(wp_width, 2.0) * RHO.get(i,i).real() + dyn_var.p_quant->add(idof, traj, 0.5 / pow(dyn_var.wp_width->get(idof, traj), 2.0) * RHO.get(i,i).real() *(dyn_var.q_aux[traj]->get(a, idof) - dyn_var.q_aux[traj]->get(i, idof))); // *(dyn_var.q->get(idof, traj) - dyn_var.q_aux[traj]->get(i, idof))); } @@ -1633,7 +1615,7 @@ void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& h *dyn_var.f += *dyn_var.f_xf; } -void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms, bp::dict params, int rotation){ +void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms, int rotation){ int itraj, i, j; int num_el = prms.num_electronic_substeps; @@ -1667,24 +1649,14 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa if(prms.assume_always_consistent){ T_new.identity(); } // Generate the half-time exponential operator - CMATRIX Hxf_old(nadi, nadi); CMATRIX Hxf(nadi, nadi); CMATRIX D(nadi, nadi); /// this is \exp[-idt/4\hbar * ( T_new.H()*Hxf(t+dt)*T_new + Hxf(t) )] - double wp_width; - if (prms.use_td_width == 1){ - wp_width = compute_wp_width(dyn_var, prms, params, dt); - } - else{ - wp_width = prms.wp_width; - } - - XF_correction(Hxf_old, dyn_var, C, wp_width, T, itraj); - //XF_correction(Hxf, dyn_var, C, prms.wp_width, T, itraj); + XF_correction(Hxf, dyn_var, C, T, itraj); + //XF_correction(Hxf, dyn_var, C, T, itraj); //Hxf = T_new.H() * Hxf * T_new; //Hxf += Hxf_old; - Hxf = Hxf_old; D = libspecialfunctions::exp_(Hxf, complex(0.0, -0.5*dt) ); if(rotation == 1){ diff --git a/src/dyn/dyn_variables.cpp b/src/dyn/dyn_variables.cpp index acf1ddbf3..cfd8ab721 100755 --- a/src/dyn/dyn_variables.cpp +++ b/src/dyn/dyn_variables.cpp @@ -203,6 +203,7 @@ void dyn_variables::allocate_shxf(){ nab_phase[itraj] = new MATRIX(nadi, ndof); } + wp_width = new MATRIX(ndof, ntraj); p_quant = new MATRIX(ndof, ntraj); VP = new MATRIX(ndof, ntraj); @@ -238,6 +239,7 @@ void dyn_variables::allocate_mqcxf(){ nab_phase[itraj] = new MATRIX(nadi, ndof); } + wp_width = new MATRIX(ndof, ntraj); p_quant = new MATRIX(ndof, ntraj); VP = new MATRIX(ndof, ntraj); f_xf = new MATRIX(ndof, ntraj); @@ -357,6 +359,7 @@ dyn_variables::dyn_variables(const dyn_variables& x){ *p_aux_old[itraj] = *x.p_aux_old[itraj]; *nab_phase[itraj] = *x.nab_phase[itraj]; } + *wp_width = *x.wp_width; *p_quant = *x.p_quant; *VP = *x.VP; @@ -373,6 +376,7 @@ dyn_variables::dyn_variables(const dyn_variables& x){ *p_aux_old[itraj] = *x.p_aux_old[itraj]; *nab_phase[itraj] = *x.nab_phase[itraj]; } + *wp_width = *x.wp_width; *p_quant = *x.p_quant; *VP = *x.VP; *f_xf = *x.f_xf; @@ -484,6 +488,7 @@ dyn_variables::~dyn_variables(){ p_aux_old.clear(); nab_phase.clear(); + delete wp_width; delete p_quant; delete VP; @@ -514,6 +519,7 @@ dyn_variables::~dyn_variables(){ p_aux_old.clear(); nab_phase.clear(); + delete wp_width; delete p_quant; delete VP; delete f_xf; diff --git a/src/dyn/dyn_variables.h b/src/dyn/dyn_variables.h index cde986815..5831c667d 100755 --- a/src/dyn/dyn_variables.h +++ b/src/dyn/dyn_variables.h @@ -388,6 +388,14 @@ class dyn_variables{ vector */ vector nab_phase; + + /** + Wave packet widths based on the Gaussian approximation + + Options: + MATRIX(ndof, ntraj) + */ + MATRIX* wp_width; /** Quantum momenta defined as (-1) * \nabla_nuc |\chi| / |\chi| @@ -441,6 +449,13 @@ class dyn_variables{ vector tcnbra_ekin; + ///================= Misc =================== + /** + The current MD time step + */ + int timestep; + + ///====================== In dyn_variables.cpp ===================== @@ -476,6 +491,7 @@ class dyn_variables{ MATRIX get_coords(){ return *q; } MATRIX get_momenta(){ return *p; } MATRIX get_forces(){ return *f; } + MATRIX get_wp_width(){ return *wp_width; } MATRIX get_p_quant(){ return *p_quant; } MATRIX get_VP(){ return *VP; } MATRIX get_f_xf(){ return *f_xf; } @@ -483,7 +499,14 @@ class dyn_variables{ MATRIX get_momenta_aux(int i){ return *p_aux[i]; } MATRIX get_nab_phase(int i){ return *nab_phase[i]; } - + void get_current_timestep(bp::dict params){ + std::string key; + for(int i=0;i(params.keys()[i]); + if(key=="timestep") { timestep = bp::extract(params.values()[i]); } + else {continue;} + } + } diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index 941deac59..8cc238e3b 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -247,6 +247,7 @@ void export_dyn_variables_objects(){ .def("get_coords", &dyn_variables::get_coords) .def("get_momenta", &dyn_variables::get_momenta) .def("get_forces", &dyn_variables::get_forces) + .def("get_wp_width", &dyn_variables::get_wp_width) .def("get_p_quant", &dyn_variables::get_p_quant) .def("get_VP", &dyn_variables::get_VP) .def("get_f_xf", &dyn_variables::get_f_xf) diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index 5d02f42bd..14eee51fd 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -305,12 +305,12 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): E_ij = <|E_i - E_j|>. It is needed when dephasing_informed option is used - * **dyn_params["wp_width"]** ( double ): A width of a Gaussian function as an approximation to adiabatic wave packets. [ default: 0.3 Bohr ] + * **dyn_params["wp_width"]** ( MATRIX(ndof, 1) ): A width of a Gaussian function as an approximation to adiabatic wave packets. This value is used for the initial conditions when the td Gaussian approximation is used, that is, use_td_width == 1 Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` - * **dyn_params["wp_v"]** ( double ): The velocity of Gaussian wave packet following potential-free td Gaussian [ default: 0.0 a.u. ] + * **dyn_params["wp_v"]** ( MATRIX(ndof,1) ): The velocity of Gaussian wave packet following potential-free td Gaussian Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` @@ -561,7 +561,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): "decoherence_rates":MATRIX(nstates, nstates), "ave_gaps":MATRIX(nstates,nstates), "schwartz_decoherence_inv_alpha": MATRIX(nstates, 1), - "wp_width":0.3, "wp_v":0, "coherence_threshold":0.01, "e_mask": 0.0001, + "wp_width":MATRIX(dyn_var.ndof, 1), "wp_v":MATRIX(dyn_var.ndof, 1), "coherence_threshold":0.01, "e_mask": 0.0001, "use_xf_force": 0, "project_out_aux": 0, "tp_algo": 1, "use_td_width": 0 } ) From e359c664ca8d5bb188551700a533cd167c824b84 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sun, 21 Jan 2024 10:03:43 -0500 Subject: [PATCH 31/47] Use the dot product between force and momentum for the collapsing --- src/dyn/dyn_decoherence_methods.cpp | 41 +++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 11 deletions(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index dd77c49e9..69de2a541 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1166,7 +1166,7 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; int a = dyn_var.act_states[traj]; - + MATRIX p_real(ndof, 1); MATRIX p_aux_temp(ndof, 1); @@ -1218,11 +1218,21 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn // Check the turning point if (is_first[i] == 0 and prms.tp_algo != 0){ - double temp = 0.0; - for(int idof=0; idofget(idof,0) ); + } + double dp_old = 0.0; + double dp_new = 0.0; + + for(int idof=0; idofget(idof,0) ); + } + double dp_old = 0.0; + double dp_new = 0.0; + for(int idof=0; idof Date: Sun, 21 Jan 2024 22:13:14 -0500 Subject: [PATCH 32/47] Add new td Gaussian width schemes --- src/dyn/Dynamics.cpp | 15 ++++++++ src/dyn/dyn_decoherence_methods.cpp | 56 +++++++++++++++++++++++++++++ src/libra_py/dynamics/tsh/save.py | 10 ++++++ 3 files changed, 81 insertions(+) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 5a867eac0..4d38da35d 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -598,6 +598,21 @@ void update_wp_width(dyn_variables& dyn_var, dyn_control_params& prms){ } } } + else if (prms.use_td_width == 2){ + double elapsed_time, s2; + + elapsed_time = prms.dt*dyn_var.timestep; + + for(int itraj=0; itrajget(idof,0), 2.0); + dyn_var.wp_width->set(idof, itraj, 4*M_PI/(s2*dyn_var.p->get(idof, itraj)) ); + } + } + } + else if (prms.use_td_width == 3){ + ;; + } } /* diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 69de2a541..f4e6dd1df 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1108,6 +1108,56 @@ for(traj = 0; traj < ntraj; traj++){ }// traj } +void td_width_aux(dyn_variables& dyn_var){ + int ntraj = dyn_var.ntraj; + int nadi = dyn_var.nadi; + int ndof = dyn_var.ndof; + + double width_temp; + + dyn_var.wp_width->set(-1, -1, 0.0); + for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; + vector& is_first = dyn_var.is_first[traj]; + + // wp_width is computed by pairwise widths based on auxiliary trajectories + MATRIX sum_inv_w2(ndof, 1); + MATRIX w2_temp(ndof, 1); + + int check_mixing = 0; + for(int i=0; i=j){continue;} + + if(is_mixed[i] == 0 or is_mixed[j] == 0){continue; } + check_mixing = 1; + + if(is_first[i] == 1 or is_first[j] == 1){ + w2_temp.set(-1, 1.0e+10); // At initial, an auxiliary pair does not contribute to wp_width + } + else{ + for(int idof=0; idofget(i, idof) - dyn_var.q_aux[traj]->get(j, idof))/ + fabs(dyn_var.p_aux[traj]->get(i, idof) - dyn_var.p_aux[traj]->get(j, idof)) ); + } + } + + for(int idof=0; idofset(idof, traj, sqrt((nadi - 1)* 1.0/sum_inv_w2.get(idof, 0)) ); + } + } + + } // traj +} + void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms){ /** \brief The generic framework of the SHXF (Surface Hopping based on eXact Factorization) method of @@ -1273,6 +1323,9 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn } }//traj + // Compute the td width based on auxiliary trajectories + if(prms.use_td_width == 3){ td_width_aux(dyn_var);} + // Propagate the spatial derivative of phases for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; @@ -1493,6 +1546,9 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy }//traj + // Compute the td width based on auxiliary trajectories + if(prms.use_td_width == 3){ td_width_aux(dyn_var);} + // Propagate the spatial derivative of phases; the E-based approximation is used for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; diff --git a/src/libra_py/dynamics/tsh/save.py b/src/libra_py/dynamics/tsh/save.py index 7868119f3..00ee7ac78 100755 --- a/src/libra_py/dynamics/tsh/save.py +++ b/src/libra_py/dynamics/tsh/save.py @@ -182,6 +182,10 @@ def init_tsh_data(saver, output_level, _nsteps, _ntraj, _ndof, _nadi, _ndia): # Trajectory-resolved diabatic TD-SE amplitudes if "Cdia" in saver.keywords: # and "Cdia" in saver.np_data.keys(): saver.add_dataset("Cdia", (_nsteps, _ntraj, _ndia), "C") + + # Trajectory-resolved quantum momenta + if "wp_width" in saver.keywords: # and "p_quant" in saver.np_data.keys(): + saver.add_dataset("wp_width", (_nsteps, _ntraj, _ndof), "R") # Trajectory-resolved quantum momenta if "p_quant" in saver.keywords: # and "p_quant" in saver.np_data.keys(): @@ -622,6 +626,12 @@ def save_hdf5_3D_new(saver, i, dyn_var, txt_type=0): if "Cdia" in saver.keywords and "Cdia" in saver.np_data.keys(): Cdia = dyn_var.get_ampl_dia() saver.save_matrix(t, "Cdia", Cdia.T()) + + # Wavepacket width + # Format: saver.add_dataset("wp_width", (_nsteps, _ntraj, _dof), "R") + if "wp_width" in saver.keywords and "wp_width" in saver.np_data.keys(): + wp_width = dyn_var.get_wp_width() + saver.save_matrix(t, "wp_width", wp_width.T()) # Trajectory-resolved quantum momenta # Format: saver.add_dataset("p_quant", (_nsteps, _ntraj, _dof), "R") From f64d3444cd4148d13b7acfc3950efd0e76f2ad17 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sun, 21 Jan 2024 22:22:05 -0500 Subject: [PATCH 33/47] Exclude projecting out the active state in MQCXF --- src/dyn/dyn_decoherence_methods.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index f4e6dd1df..8cf785a09 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1403,6 +1403,8 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy vector& is_keep = dyn_var.is_keep[traj]; MATRIX& p_aux = *dyn_var.p_aux[traj]; MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; + + int a = dyn_var.act_states[traj]; MATRIX p_real(ndof, 1); MATRIX p_aux_temp(ndof, 1); @@ -1440,7 +1442,7 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy } if (alpha < 0.0){ alpha = 0.0; - if (prms.project_out_aux == 1){ + if (prms.project_out_aux == 1 and i!=a){ project_out(*dyn_var.ampl_adi, traj, i); xf_init_AT(dyn_var, traj, -1); cout << "Project out a classically forbidden state " << i << " on traj " << traj < Date: Fri, 2 Feb 2024 02:01:11 -0500 Subject: [PATCH 34/47] Full propagation before reordering --- src/dyn/Dynamics.cpp | 6 +++--- src/dyn/dyn_decoherence_methods.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 4d38da35d..68097c459 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1272,9 +1272,9 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, propagate_electronic(dyn_var, ham, ham_aux, prms); - if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ - propagate_half_xf(dyn_var, ham, prms, 1); - } + //if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ + // propagate_half_xf(dyn_var, ham, prms, 1); + //} } diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 8cf785a09..7eb5acb47 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1735,7 +1735,7 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa //Hxf = T_new.H() * Hxf * T_new; //Hxf += Hxf_old; - D = libspecialfunctions::exp_(Hxf, complex(0.0, -0.5*dt) ); + D = libspecialfunctions::exp_(Hxf, complex(0.0, -dt) ); if(rotation == 1){ D = T_new * D * T_new.H(); } From f16d6bc6acbfaef9065fc5c1c5bf2535a4d469a3 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Fri, 2 Feb 2024 02:55:34 -0500 Subject: [PATCH 35/47] Adjust the nabla phase according to the reordering --- src/dyn/Dynamics.cpp | 5 +- src/dyn/dyn_decoherence.h | 1 + src/dyn/dyn_decoherence_methods.cpp | 146 +++++++++++++++++++++------- src/dyn/dyn_variables.cpp | 10 ++ src/dyn/dyn_variables.h | 8 ++ 5 files changed, 134 insertions(+), 36 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 68097c459..ba412f8a5 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1276,7 +1276,10 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, // propagate_half_xf(dyn_var, ham, prms, 1); //} } - + + if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ + update_nab_phase(dyn_var, ham, prms); + } // Recompute density matrices in response to the updated amplitudes dyn_var.update_amplitudes(prms, ham); diff --git a/src/dyn/dyn_decoherence.h b/src/dyn/dyn_decoherence.h index 35b9bebe3..9ae19c2d1 100755 --- a/src/dyn/dyn_decoherence.h +++ b/src/dyn/dyn_decoherence.h @@ -80,6 +80,7 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms); // For MQCXF // XF propagation +void update_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms); void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev); void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms, int rotation); void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj); diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 7eb5acb47..a0d0f3848 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1326,33 +1326,34 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn // Compute the td width based on auxiliary trajectories if(prms.use_td_width == 3){ td_width_aux(dyn_var);} - // Propagate the spatial derivative of phases - for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; - vector& is_first = dyn_var.is_first[traj]; - MATRIX& p_aux = *dyn_var.p_aux[traj]; - MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; - MATRIX& nab_phase = *dyn_var.nab_phase[traj]; - - for(int i=0; i& is_mixed = dyn_var.is_mixed[traj]; + // vector& is_first = dyn_var.is_first[traj]; + // MATRIX& p_aux = *dyn_var.p_aux[traj]; + // MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; + // MATRIX& nab_phase = *dyn_var.nab_phase[traj]; + + // for(int i=0; iget(0,0) << " " << dyn_var.q_aux[0]->get(0,0) << " " << dyn_var.q_aux[0]->get(1,0) // << " " << dyn_var.p_aux[0]->get(0,0) << " " << dyn_var.p_aux[0]->get(1,0) <& is_mixed = dyn_var.is_mixed[traj]; + // vector& is_first = dyn_var.is_first[traj]; + // MATRIX& nab_phase = *dyn_var.nab_phase[traj]; + // + // CMATRIX E(nadi, nadi); + // E = ham.children[traj]->get_ham_adi(); + + // MATRIX p_real(ndof, 1); + // p_real = dyn_var.p->col(traj); + // double Ekin = compute_kinetic_energy(p_real, invM); + + // for(int i=0; iget(idof,traj)/(Ekin + prms.e_mask)); + // }//idof + // } + // }//i + //} // traj +} + +void update_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms){ + int ntraj = dyn_var.ntraj; + int nadi = dyn_var.nadi; + int ndof = dyn_var.ndof; + for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; vector& is_first = dyn_var.is_first[traj]; + MATRIX& p_aux = *dyn_var.p_aux[traj]; + MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; MATRIX& nab_phase = *dyn_var.nab_phase[traj]; - - CMATRIX E(nadi, nadi); - E = ham.children[traj]->get_ham_adi(); + MATRIX& nab_phase_old = *dyn_var.nab_phase_old[traj]; + + CMATRIX T_new(*dyn_var.proj_adi[traj]); + if(prms.assume_always_consistent){ T_new.identity(); } + + vector P_AUX; vector P_AUX_OLD; + for(int idof=0; idofcol(traj); - double Ekin = compute_kinetic_energy(p_real, invM); + if (prms.decoherence_algo==5){ + for(int i=0; iget(idof,traj)/(Ekin + prms.e_mask)); - }//idof - } - }//i + for(int idof=0; idofget_ham_adi(); + E = T_new.H() * E * T_new; + + MATRIX p_real(ndof, 1); + p_real = dyn_var.p->col(traj); + double Ekin = compute_kinetic_energy(p_real, invM); + + for(int i=0; iget(idof,traj)/(Ekin + prms.e_mask)); + }//idof + } + }//i + + } } // traj + } void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, CMATRIX& T, int traj){ diff --git a/src/dyn/dyn_variables.cpp b/src/dyn/dyn_variables.cpp index cfd8ab721..9c1e91a6f 100755 --- a/src/dyn/dyn_variables.cpp +++ b/src/dyn/dyn_variables.cpp @@ -195,12 +195,14 @@ void dyn_variables::allocate_shxf(){ p_aux = vector(ntraj); p_aux_old = vector(ntraj); nab_phase = vector(ntraj); + nab_phase_old = vector(ntraj); for(int itraj=0; itraj(ntraj); p_aux_old = vector(ntraj); nab_phase = vector(ntraj); + nab_phase_old = vector(ntraj); for(int itraj=0; itraj nab_phase; + /** + Spatial derivative of the phase of coefficients of state-wise auxiliary trajectories + + Options: + vector + */ + vector nab_phase_old; + /** Wave packet widths based on the Gaussian approximation From 36cec5bb3c8895032e2a4237e8e6cb5f3f951652 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Sun, 4 Feb 2024 21:06:23 -0500 Subject: [PATCH 36/47] Apply the half propagations and reordering --- src/dyn/Dynamics.cpp | 12 ++++-------- src/dyn/dyn_decoherence_methods.cpp | 2 +- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index ba412f8a5..e023f66cc 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1272,19 +1272,15 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, propagate_electronic(dyn_var, ham, ham_aux, prms); - //if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ - // propagate_half_xf(dyn_var, ham, prms, 1); - //} - } - - if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ - update_nab_phase(dyn_var, ham, prms); + if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ + update_nab_phase(dyn_var, ham, prms); + propagate_half_xf(dyn_var, ham, prms, 0); + } } // Recompute density matrices in response to the updated amplitudes dyn_var.update_amplitudes(prms, ham); dyn_var.update_density_matrix(prms, ham, 1); - vector old_states( dyn_var.act_states); diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index a0d0f3848..6aaa8a5fa 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1811,7 +1811,7 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa //Hxf = T_new.H() * Hxf * T_new; //Hxf += Hxf_old; - D = libspecialfunctions::exp_(Hxf, complex(0.0, -dt) ); + D = libspecialfunctions::exp_(Hxf, complex(0.0, -0.5*dt) ); if(rotation == 1){ D = T_new * D * T_new.H(); } From a200cf6cd0a1bf2f60b1623fca809fd82ee6150c Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Tue, 6 Feb 2024 01:22:50 -0500 Subject: [PATCH 37/47] removed Gaussian_methods from the libra_py/__init__.py file --- src/libra_py/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/libra_py/__init__.py b/src/libra_py/__init__.py index 7316cf3f3..0f419804c 100755 --- a/src/libra_py/__init__.py +++ b/src/libra_py/__init__.py @@ -28,7 +28,6 @@ "fix_motion", "ft", "gaussian_kernel_algorithm", - "Gaussian_methods", "hpc_utils", "hungarian", "init_ensembles", From af69fe6550524167d06200de2414dd73f9be05d2 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Thu, 8 Feb 2024 05:01:27 -0500 Subject: [PATCH 38/47] Remove the thermostating stuff --- src/dyn/Dynamics.cpp | 30 ++---------------------------- 1 file changed, 2 insertions(+), 28 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index e023f66cc..620a5c668 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1403,38 +1403,12 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, } // SHXF else if(prms.decoherence_algo==5){ - if(prms.rep_tdse==1){ - shxf(dyn_var, ham, ham_aux, prms); - - if(prms.ensemble==1){ - for(idof=0; idofscale(i, dof, therm[traj].vel_scale(1.0*prms.dt)); - }// i - }// traj - }//idof - } - } + if(prms.rep_tdse==1){shxf(dyn_var, ham, ham_aux, prms);} else{ cout<<"ERROR: SHXF requires rep_tdse = 1\nExiting now...\n"; exit(0); } } // MQCXF else if(prms.decoherence_algo==6){ - if(prms.rep_tdse==1){ - mqcxf(dyn_var, ham, ham_aux, prms); - - if(prms.ensemble==1){ - for(idof=0; idofscale(i, dof, therm[traj].vel_scale(1.0*prms.dt)); - }// i - }// traj - }//idof - } - } + if(prms.rep_tdse==1){mqcxf(dyn_var, ham, ham_aux, prms);} else{ cout<<"ERROR: MQCXF requires rep_tdse = 1\nExiting now...\n"; exit(0); } } From 1018adf9b628413dbbb5ac91a062e90754141828 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Thu, 8 Feb 2024 05:23:03 -0500 Subject: [PATCH 39/47] Remove the redundant calls of T in the XF propagation --- src/dyn/Dynamics.cpp | 4 ++-- src/dyn/dyn_decoherence.h | 2 +- src/dyn/dyn_decoherence_methods.cpp | 24 +++++------------------- 3 files changed, 8 insertions(+), 22 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 620a5c668..2a05e0a54 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1267,14 +1267,14 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, // basis re-projection matrices for(i=0; i(0.0, dyn_var.nab_phase[traj]->get(i, idof)) ); } } - F = T * F * T.H(); Ham += -invM.get(idof,0) * dyn_var.p_quant->get(idof, traj)*(RHO * F - F * RHO); } @@ -1768,7 +1767,7 @@ void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& h *dyn_var.f += *dyn_var.f_xf; } -void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms, int rotation){ +void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms){ int itraj, i, j; int num_el = prms.num_electronic_substeps; @@ -1797,29 +1796,16 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa nHamiltonian* ham = Ham.children[traj1]; - CMATRIX T(*dyn_var.proj_adi[itraj]); T.load_identity(); - CMATRIX T_new(*dyn_var.proj_adi[itraj]); - if(prms.assume_always_consistent){ T_new.identity(); } - // Generate the half-time exponential operator CMATRIX Hxf(nadi, nadi); - CMATRIX D(nadi, nadi); /// this is \exp[-idt/4\hbar * ( T_new.H()*Hxf(t+dt)*T_new + Hxf(t) )] - - XF_correction(Hxf, dyn_var, C, T, itraj); - //XF_correction(Hxf, dyn_var, C, T, itraj); + CMATRIX D(nadi, nadi); // this is \exp[ -idt/2\hbar * Hxf( C(t) ) ] - //Hxf = T_new.H() * Hxf * T_new; - //Hxf += Hxf_old; + XF_correction(Hxf, dyn_var, C, itraj); D = libspecialfunctions::exp_(Hxf, complex(0.0, -0.5*dt) ); - if(rotation == 1){ - D = T_new * D * T_new.H(); - } C = D * C; -// *dyn_var.proj_adi[itraj] = T_new; - // Insert the propagated result back for(int st=0; st Date: Thu, 8 Feb 2024 07:06:31 -0500 Subject: [PATCH 40/47] Separate the XF Hamiltonian calculation from the propagation --- src/dyn/Dynamics.cpp | 2 + src/dyn/dyn_decoherence.h | 1 + src/dyn/dyn_decoherence_methods.cpp | 91 +++++++++++++++-------------- src/dyn/dyn_variables.cpp | 10 ++++ src/dyn/dyn_variables.h | 8 +++ src/dyn/libdyn.cpp | 4 ++ 6 files changed, 72 insertions(+), 44 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 2a05e0a54..5eada493e 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1267,6 +1267,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, // basis re-projection matrices for(i=0; i& // Independent-trajectory XF methods void xf_hop_reset(dyn_variables& dyn_var, vector& accepted_states, vector& initial_states); +void update_ham_xf(dyn_variables& dyn_var); void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms); // For SHXF void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms); // For MQCXF diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index bf7c62075..7f6431f10 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1650,45 +1650,58 @@ void update_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_par } -void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, int traj){ +void update_ham_xf(dyn_variables& dyn_var){ int ndof = dyn_var.ndof; + int ntraj = dyn_var.ntraj; int nst = dyn_var.nadi; MATRIX& invM = *dyn_var.iM; - - vector& is_mixed = dyn_var.is_mixed[traj]; // Construct and transform the density matrix + CMATRIX C(nst, 1); + CMATRIX Coeff(nst, ntraj); CMATRIX RHO(nst, nst); - RHO = C * C.H(); + + Coeff = *dyn_var.ampl_adi; - // Compute quantum momenta - int a = dyn_var.act_states[traj]; + // temp for the XF Hamiltonian + CMATRIX iphi(nst, nst); + CMATRIX Ham_xf(nst, nst); - dyn_var.p_quant->set(-1, traj, 0.0); - for(int i=0; iadd(idof, traj, 0.5 / pow(dyn_var.wp_width->get(idof, traj), 2.0) * RHO.get(i,i).real() - *(dyn_var.q_aux[traj]->get(a, idof) - dyn_var.q_aux[traj]->get(i, idof))); - // *(dyn_var.q->get(idof, traj) - dyn_var.q_aux[traj]->get(i, idof))); - } - } - } + for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; - // Add the XF-based decoherence correction - for(int idof=0; idofset(-1, traj, 0.0); for(int i=0; i(0.0, dyn_var.nab_phase[traj]->get(i, idof)) ); + for(int idof=0; idofadd(idof, traj, 0.5 / pow(dyn_var.wp_width->get(idof, traj), 2.0) * RHO.get(i,i).real() + *(dyn_var.q_aux[traj]->get(a, idof) - dyn_var.q_aux[traj]->get(i, idof))); + // *(dyn_var.q->get(idof, traj) - dyn_var.q_aux[traj]->get(i, idof))); + } } } - Ham += -invM.get(idof,0) * dyn_var.p_quant->get(idof, traj)*(RHO * F - F * RHO); - } + // Add the XF-based decoherence correction + Ham_xf.set(-1, -1, complex(0.0, 0.0)); + for(int idof=0; idof(0.0, dyn_var.nab_phase[traj]->get(i, idof)) ); + } + } + + Ham_xf += -invM.get(idof,0) * dyn_var.p_quant->get(idof, traj)*(RHO * iphi - iphi * RHO); + } + *dyn_var.ham_xf[traj] = Ham_xf; + } //traj } void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev){ @@ -1705,9 +1718,9 @@ void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& h CMATRIX Coeff(nst, ntraj); // termporaries for nabla_phase and adiabatic force - vector F; + vector phi; for(int idof=0; idofcol(traj); double Ekin = compute_kinetic_energy(p_real, invM); - // Compute F for each dof + // Compute phi for each dof vector& is_mixed = dyn_var.is_mixed[traj]; for(int idof=0; idof (0.0, 0.0)); + phi[idof].set(-1,-1, complex (0.0, 0.0)); for(int i=0; i (dyn_var.nab_phase[traj]->get(i, idof), 0.0) ); + phi[idof].set(i,i,complex (dyn_var.nab_phase[traj]->get(i, idof), 0.0) ); } } } // Save vector potential (contribution from NAC is neglected) for(int idof=0; idofset(idof, traj, temp.get(0,0).real()); } // Original form of the decoherence force for(int idof=0; idofadd(idof, traj, -2.0*invM.get(jdof,0)*dyn_var.p_quant->get(jdof, traj)* (dyn_var.VP->get(jdof, traj)*dyn_var.VP->get(idof, traj) - temp.get(0,0).real() ) ); } } - //// Energy-conserving treatment - //for(int idof=0; idofget(jdof, traj)* - // ( dyn_var.VP->get(jdof, traj)*Epot - (C.H()*F[jdof]*E*C).get(0,0).real() ); - // } - // dyn_var.f_xf->set(idof, traj, temp / Ekin * dyn_var.p->get(idof, traj) ); - //} } //traj // Add the XF contribution @@ -1787,6 +1791,9 @@ void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_pa CMATRIX Coeff(nst, ntraj); Coeff = *dyn_var.ampl_adi; + + CMATRIX Hxf(nadi, nadi); + CMATRIX D(nadi, nadi); // this is \exp[ -idt/2\hbar * Hxf( C(t) ) ] for(itraj=0; itraj(0.0, -0.5*dt) ); C = D * C; diff --git a/src/dyn/dyn_variables.cpp b/src/dyn/dyn_variables.cpp index 9c1e91a6f..e65cb89bd 100755 --- a/src/dyn/dyn_variables.cpp +++ b/src/dyn/dyn_variables.cpp @@ -196,6 +196,7 @@ void dyn_variables::allocate_shxf(){ p_aux_old = vector(ntraj); nab_phase = vector(ntraj); nab_phase_old = vector(ntraj); + ham_xf = vector(ntraj); for(int itraj=0; itraj(ntraj); nab_phase = vector(ntraj); nab_phase_old = vector(ntraj); + ham_xf = vector(ntraj); for(int itraj=0; itraj nab_phase_old; + /** + XF Hamiltonian + + Options: + vector + */ + vector ham_xf; + /** Wave packet widths based on the Gaussian approximation diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index 8cc238e3b..af74caf00 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -366,6 +366,10 @@ void export_dyn_decoherence_objects(){ void (*expt_xf_hop_reset) (dyn_variables& dyn_var, vector& accepted_states, vector& initial_states) = &xf_hop_reset; def("xf_hop_reset", expt_xf_hop_reset); + + void (*expt_update_ham_xf) + (dyn_variables& dyn_var) = &update_ham_xf; + def("update_ham_xf", expt_update_ham_xf); void (*expt_shxf_v1) (dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms) = &shxf; From ecaa49bc85ce4fb9554da52e96fa15469ef2f5fc Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Thu, 8 Feb 2024 08:05:38 -0500 Subject: [PATCH 41/47] Compute the phase gradient in the usual places and only rotate it before the propagation --- src/dyn/Dynamics.cpp | 2 +- src/dyn/dyn_decoherence.h | 2 +- src/dyn/dyn_decoherence_methods.cpp | 165 +++++++++++----------------- src/dyn/dyn_variables.cpp | 10 +- 4 files changed, 74 insertions(+), 105 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 5eada493e..45566e26f 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -1274,7 +1274,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, propagate_electronic(dyn_var, ham, ham_aux, prms); if(prms.decoherence_algo == 5 or prms.decoherence_algo == 6){ - update_nab_phase(dyn_var, ham, prms); + rotate_nab_phase(dyn_var, ham, prms); update_ham_xf(dyn_var); propagate_half_xf(dyn_var, ham, prms); } diff --git a/src/dyn/dyn_decoherence.h b/src/dyn/dyn_decoherence.h index e1c53a4e0..164332ca6 100755 --- a/src/dyn/dyn_decoherence.h +++ b/src/dyn/dyn_decoherence.h @@ -81,7 +81,7 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms); // For MQCXF // XF propagation -void update_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms); +void rotate_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms); void update_forces_xf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev); void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms); void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj); diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 7f6431f10..5879c0a46 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1326,27 +1326,28 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn // Compute the td width based on auxiliary trajectories if(prms.use_td_width == 3){ td_width_aux(dyn_var);} - //// Propagate the spatial derivative of phases - //for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; - // vector& is_first = dyn_var.is_first[traj]; - // MATRIX& p_aux = *dyn_var.p_aux[traj]; - // MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; - // MATRIX& nab_phase = *dyn_var.nab_phase[traj]; - - // for(int i=0; i& is_mixed = dyn_var.is_mixed[traj]; + vector& is_first = dyn_var.is_first[traj]; + MATRIX& p_aux = *dyn_var.p_aux[traj]; + MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; + MATRIX& nab_phase = *dyn_var.nab_phase[traj]; + MATRIX& nab_phase_old = *dyn_var.nab_phase_old[traj]; + + for(int i=0; iget(0,0) << " " << dyn_var.q_aux[0]->get(0,0) << " " << dyn_var.q_aux[0]->get(1,0) // << " " << dyn_var.p_aux[0]->get(0,0) << " " << dyn_var.p_aux[0]->get(1,0) <& is_mixed = dyn_var.is_mixed[traj]; - // vector& is_first = dyn_var.is_first[traj]; - // MATRIX& nab_phase = *dyn_var.nab_phase[traj]; - // - // CMATRIX E(nadi, nadi); - // E = ham.children[traj]->get_ham_adi(); - - // MATRIX p_real(ndof, 1); - // p_real = dyn_var.p->col(traj); - // double Ekin = compute_kinetic_energy(p_real, invM); - - // for(int i=0; iget(idof,traj)/(Ekin + prms.e_mask)); - // }//idof - // } - // }//i - //} // traj + // Propagate the spatial derivative of phases; the E-based approximation is used + for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; + vector& is_first = dyn_var.is_first[traj]; + MATRIX& nab_phase = *dyn_var.nab_phase[traj]; + + CMATRIX E(nadi, nadi); + E = ham.children[traj]->get_ham_adi(); + + MATRIX p_real(ndof, 1); + p_real = dyn_var.p->col(traj); + double Ekin = compute_kinetic_energy(p_real, invM); + + for(int i=0; iget(idof,traj)/(Ekin + prms.e_mask)); + }//idof + } + }//i + } // traj } -void update_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms){ +void rotate_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms){ int ntraj = dyn_var.ntraj; - int nadi = dyn_var.nadi; + int nst = dyn_var.nadi; int ndof = dyn_var.ndof; + CMATRIX phi(nst, nst); + for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; - vector& is_first = dyn_var.is_first[traj]; - MATRIX& p_aux = *dyn_var.p_aux[traj]; - MATRIX& p_aux_old = *dyn_var.p_aux_old[traj]; - MATRIX& nab_phase = *dyn_var.nab_phase[traj]; - MATRIX& nab_phase_old = *dyn_var.nab_phase_old[traj]; - + CMATRIX T_new(*dyn_var.proj_adi[traj]); if(prms.assume_always_consistent){ T_new.identity(); } - - vector P_AUX; vector P_AUX_OLD; + + // Set a diagonal matrix of nabla_phase for each dof for(int idof=0; idof(dyn_var.nab_phase[traj]->get(i, idof), 0.0) ); } - }//i - } - else{ - MATRIX& invM = *dyn_var.iM; - CMATRIX E(nadi, nadi); - E = ham.children[traj]->get_ham_adi(); - E = T_new.H() * E * T_new; - - MATRIX p_real(ndof, 1); - p_real = dyn_var.p->col(traj); - double Ekin = compute_kinetic_energy(p_real, invM); + } + phi = T_new.H() * phi * T_new; + + for(int i=0; iset(i, idof, phi.get(i,i).real() ); + } + } + }//idof - for(int i=0; iget(idof,traj)/(Ekin + prms.e_mask)); - }//idof + dyn_var.nab_phase_old[traj]->set(i, idof, dyn_var.nab_phase[traj]->get(i, idof) ); + } } - }//i - + } } + } // traj } diff --git a/src/dyn/dyn_variables.cpp b/src/dyn/dyn_variables.cpp index e65cb89bd..38bfffbae 100755 --- a/src/dyn/dyn_variables.cpp +++ b/src/dyn/dyn_variables.cpp @@ -235,7 +235,7 @@ void dyn_variables::allocate_mqcxf(){ p_aux = vector(ntraj); p_aux_old = vector(ntraj); nab_phase = vector(ntraj); - nab_phase_old = vector(ntraj); + //nab_phase_old = vector(ntraj); ham_xf = vector(ntraj); for(int itraj=0; itraj Date: Fri, 9 Feb 2024 01:05:27 -0500 Subject: [PATCH 42/47] Move the Subotnik td width calculation to update_wp_width --- src/dyn/Dynamics.cpp | 47 +++++++++++++++++++++++- src/dyn/dyn_decoherence_methods.cpp | 56 ----------------------------- 2 files changed, 46 insertions(+), 57 deletions(-) diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 45566e26f..5df9a3b20 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -611,7 +611,52 @@ void update_wp_width(dyn_variables& dyn_var, dyn_control_params& prms){ } } else if (prms.use_td_width == 3){ - ;; + int nadi = dyn_var.nadi; + + // wp_width is computed by pairwise widths based on auxiliary trajectories + dyn_var.wp_width->set(-1, -1, 0.0); + MATRIX sum_inv_w2(ndof, 1); + MATRIX w2_temp(ndof, 1); + + for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; + vector& is_first = dyn_var.is_first[traj]; + + sum_inv_w2.set(-1, 0, 0.0); + w2_temp.set(-1, 0, 0.0); + + int check_mixing = 0; + for(int i=0; i=j){continue;} + + if(is_mixed[i] == 0 or is_mixed[j] == 0){continue; } + check_mixing = 1; + + if(is_first[i] == 1 or is_first[j] == 1){ + w2_temp.set(-1, 0, 1.0e+10); // At initial, an auxiliary pair does not contribute to wp_width + } + else{ + for(int idof=0; idofget(i, idof) - dyn_var.q_aux[traj]->get(j, idof))/ + fabs(dyn_var.p_aux[traj]->get(i, idof) - dyn_var.p_aux[traj]->get(j, idof)) ); + } + } + + for(int idof=0; idofset(idof, traj, sqrt((nadi - 1)* 1.0/sum_inv_w2.get(idof, 0)) ); + } + } + + } // traj } } diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index 5879c0a46..e24cca69c 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1108,56 +1108,6 @@ for(traj = 0; traj < ntraj; traj++){ }// traj } -void td_width_aux(dyn_variables& dyn_var){ - int ntraj = dyn_var.ntraj; - int nadi = dyn_var.nadi; - int ndof = dyn_var.ndof; - - double width_temp; - - dyn_var.wp_width->set(-1, -1, 0.0); - for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; - vector& is_first = dyn_var.is_first[traj]; - - // wp_width is computed by pairwise widths based on auxiliary trajectories - MATRIX sum_inv_w2(ndof, 1); - MATRIX w2_temp(ndof, 1); - - int check_mixing = 0; - for(int i=0; i=j){continue;} - - if(is_mixed[i] == 0 or is_mixed[j] == 0){continue; } - check_mixing = 1; - - if(is_first[i] == 1 or is_first[j] == 1){ - w2_temp.set(-1, 1.0e+10); // At initial, an auxiliary pair does not contribute to wp_width - } - else{ - for(int idof=0; idofget(i, idof) - dyn_var.q_aux[traj]->get(j, idof))/ - fabs(dyn_var.p_aux[traj]->get(i, idof) - dyn_var.p_aux[traj]->get(j, idof)) ); - } - } - - for(int idof=0; idofset(idof, traj, sqrt((nadi - 1)* 1.0/sum_inv_w2.get(idof, 0)) ); - } - } - - } // traj -} - void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms){ /** \brief The generic framework of the SHXF (Surface Hopping based on eXact Factorization) method of @@ -1323,9 +1273,6 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn } }//traj - // Compute the td width based on auxiliary trajectories - if(prms.use_td_width == 3){ td_width_aux(dyn_var);} - // Propagate the spatial derivative of phases for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; @@ -1550,9 +1497,6 @@ void mqcxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dy }//traj - // Compute the td width based on auxiliary trajectories - if(prms.use_td_width == 3){ td_width_aux(dyn_var);} - // Propagate the spatial derivative of phases; the E-based approximation is used for(int traj=0; traj& is_mixed = dyn_var.is_mixed[traj]; From c7dba0e8b05a624a29ac0e85b6c3bec9c56599c2 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Wed, 14 Feb 2024 01:01:35 -0500 Subject: [PATCH 43/47] Add missing parts in model definition --- src/libra_py/models/Subotnik.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/libra_py/models/Subotnik.py b/src/libra_py/models/Subotnik.py index f1143e114..7302221dd 100755 --- a/src/libra_py/models/Subotnik.py +++ b/src/libra_py/models/Subotnik.py @@ -99,7 +99,7 @@ def dumbbell_geometry(q, params, full_id=None): H01 = B * ( e_minus + (2.0 - e_plus) ) dH01 = B * C * ( e_minus - e_plus) - elif x>-Z and x=-Z and x<=Z: e_plus = math.exp(-C*(x+Z)) e_minus = math.exp(C*(x-Z)) @@ -210,7 +210,7 @@ def double_arch_geometry(q, params, full_id=None): H01 = B * ( -e_minus + e_plus) dH01 = B * C * ( -e_minus + e_plus) - elif x>-Z and x=-Z and x<=Z: e_plus = math.exp(-C*(x+Z)) e_minus = math.exp(C*(x-Z)) From 5977b62766a6363d5b6b08aee43f685dd5cd99be Mon Sep 17 00:00:00 2001 From: Alexey Akimov Date: Fri, 16 Feb 2024 21:40:00 -0500 Subject: [PATCH 44/47] Added the option to compute pDOS for spin-polarized calculations with cp2k output --- src/libra_py/packages/cp2k/methods.py | 151 ++++++++++++++++++++------ 1 file changed, 118 insertions(+), 33 deletions(-) diff --git a/src/libra_py/packages/cp2k/methods.py b/src/libra_py/packages/cp2k/methods.py index e1c59a68a..a91a6d2ed 100755 --- a/src/libra_py/packages/cp2k/methods.py +++ b/src/libra_py/packages/cp2k/methods.py @@ -1689,6 +1689,49 @@ def gaussian_function_vector(a_vec, mu_vec, sigma, num_points, x_min, x_max): return x, sum_vec + +def aux_pdos(c1, atoms, pdos_files, margin, homo_occ, orbitals_cols, sigma, npoints, ave_pdos_convolved_all, labels): + """ + c1 - index of the atom kinds + atoms - the + pdos_files (list of strings): names of the files containing pdos data + shift (double): the amount of "padding" around the minimal and maximal values of the energy scale [units: eV] + homo_occ (double): occupancy of the HOMO: 2 - for non-polarized, 1 - for polarized + orbital_cols (list of lists of integers): integers that define which columns correspond to which types of orbitals + sigma (double): broadening of the pDOS + npoints (int): how many points to use to compte the convolved pDOS + + """ + + pdos_ave = np.zeros(np.loadtxt(pdos_files[0]).shape) + + cnt = len(pdos_files) + for c2, pdos_file in enumerate(pdos_files): + pdos_mat = np.loadtxt(pdos_file) + if c2==0: + pdos_ave = np.zeros(pdos_mat.shape) + pdos_ave += pdos_mat + pdos_ave /= cnt + + pdos_ave[:,1] *= units.au2ev + e_min = np.min(pdos_ave[:,1]) - margin + e_max = np.max(pdos_ave[:,1]) + margin + homo_level = np.max(np.where(pdos_ave[:,2]==homo_occ)) + homo_energy = pdos_ave[:,1][homo_level] + + for c3, orbital_cols in enumerate(orbitals_cols): + try: + sum_pdos_ave = np.sum(pdos_ave[:,orbital_cols],axis=1) + ave_energy_grid, ave_pdos_convolved = gaussian_function_vector(sum_pdos_ave, pdos_ave[:,1], sigma,npoints, e_min, e_max) + ave_pdos_convolved_all.append(ave_pdos_convolved) + pdos_label = F"{atoms[1][c1]}, {orbitals[c3]}" + labels.append(pdos_label) + except: + pass + + return ave_energy_grid, homo_energy, ave_pdos_convolved_all + + def pdos(params): """ This function computes the weighted pdos for a set of CP2K pdos files. @@ -1705,6 +1748,7 @@ def pdos(params): * npoints (integer): The number of grid points for convolution with Gaussian functions. * sigma (float): The standard deviation value * shift (float): The amount of shifting the grid from the minimum and maximum energy values. + * is_spin_polarized (int): 0 - non-spin-polarized, 1 - spin-polarized Returns: ave_energy_grid (nparray): The average energy grid for all pdos files homo_energy (float): The value of HOMO energy in eV @@ -1716,7 +1760,10 @@ def pdos(params): # Critical parameters critical_params = [ "path_to_all_pdos", "atoms"] # Default parameters - default_params = { "orbitals_cols": [[3], range(4,7), range(7,12), range(12,19)], "orbitals": ['s', 'p', 'd', 'f'], "sigma": 0.05, "shift": 2.0, "npoints": 2000} + default_params = { "orbitals_cols": [[3], range(4,7), range(7,12), range(12,19)], + "orbitals": ['s', 'p', 'd', 'f'], + "sigma": 0.05, "shift": 2.0, "npoints": 2000, "is_spin_polarized": 0 + } # Check input comn.check_input(params, default_params, critical_params) @@ -1727,41 +1774,79 @@ def pdos(params): npoints = params['npoints'] sigma = params['sigma'] shift = params['shift'] - - labels = [] - ave_pdos_convolved_all = [] - - for c1,i in enumerate(atoms[0]): - # Finding all the pdos files - pdos_files = glob.glob(path_to_all_pdos+F'/*k{i}*.pdos') - pdos_ave = np.zeros(np.loadtxt(pdos_files[0]).shape) - cnt = len(pdos_files) - for c2, pdos_file in enumerate(pdos_files): - pdos_mat = np.loadtxt(pdos_file) - if c2==0: - pdos_ave = np.zeros(pdos_mat.shape) - pdos_ave += pdos_mat - pdos_ave /= cnt - pdos_ave[:,1] *= units.au2ev - e_min = np.min(pdos_ave[:,1])-shift - e_max = np.max(pdos_ave[:,1])+shift - homo_level = np.max(np.where(pdos_ave[:,2]==2.0)) - homo_energy = pdos_ave[:,1][homo_level] - for c3, orbital_cols in enumerate(orbitals_cols): - try: - sum_pdos_ave = np.sum(pdos_ave[:,orbital_cols],axis=1) - ave_energy_grid, ave_pdos_convolved = gaussian_function_vector(sum_pdos_ave, pdos_ave[:,1], sigma, + is_spin_polarized = params["is_spin_polarized"] + + + if is_spin_polarized == 0: # non-polarized case + homo_occ = 2.0 + labels = [] + ave_pdos_convolved_all = [] + + for c1,i in enumerate(atoms[0]): + # Finding all the pdos files + pdos_files = glob.glob(path_to_all_pdos+F'/*k{i}*.pdos') + + ave_energy_grid, homo_energy, ave_pdos_convolved_all = aux_pdos(c1, atoms, pdos_files, shift, homo_occ, orbitals_cols, sigma, npoints, ave_pdos_convolved_all, labels) + + """ + # Finding all the pdos files + pdos_files = glob.glob(path_to_all_pdos+F'/*k{i}*.pdos') + pdos_ave = np.zeros(np.loadtxt(pdos_files[0]).shape) + + cnt = len(pdos_files) + for c2, pdos_file in enumerate(pdos_files): + pdos_mat = np.loadtxt(pdos_file) + if c2==0: + pdos_ave = np.zeros(pdos_mat.shape) + pdos_ave += pdos_mat + pdos_ave /= cnt + pdos_ave[:,1] *= units.au2ev + e_min = np.min(pdos_ave[:,1])-shift + e_max = np.max(pdos_ave[:,1])+shift + homo_level = np.max(np.where(pdos_ave[:,2]==2.0)) + homo_energy = pdos_ave[:,1][homo_level] + for c3, orbital_cols in enumerate(orbitals_cols): + try: + sum_pdos_ave = np.sum(pdos_ave[:,orbital_cols],axis=1) + ave_energy_grid, ave_pdos_convolved = gaussian_function_vector(sum_pdos_ave, pdos_ave[:,1], sigma, npoints, e_min, e_max) - ave_pdos_convolved_all.append(ave_pdos_convolved) - pdos_label = atoms[1][c1]+F', {orbitals[c3]}' - labels.append(pdos_label) - except: - pass + ave_pdos_convolved_all.append(ave_pdos_convolved) + pdos_label = atoms[1][c1]+F', {orbitals[c3]}' + labels.append(pdos_label) + except: + pass + """ + + ave_pdos_convolved_all = np.array(ave_pdos_convolved_all) + ave_pdos_convolved_total = np.sum(ave_pdos_convolved_all, axis=0) + + return ave_energy_grid, homo_energy, ave_pdos_convolved_all, labels, ave_pdos_convolved_total + + else: # spin-polarized case + homo_occ = 1.0 + + labels_alp, labels_bet = [], [] + ave_pdos_convolved_all_alp, ave_pdos_convolved_all_bet = [], [] + + for c1,i in enumerate(atoms[0]): + # Finding all the pdos files + pdos_files_alp = glob.glob(path_to_all_pdos+F'/*ALPHA*k{i}*.pdos') + pdos_files_bet = glob.glob(path_to_all_pdos+F'/*BETA*k{i}*.pdos') + + ave_energy_grid_alp, homo_energy_alp, ave_pdos_convolved_all_alp = aux_pdos(c1, atoms, pdos_files_alp, shift, homo_occ, orbitals_cols, sigma, npoints, ave_pdos_convolved_all_alp, labels_alp) + + ave_energy_grid_bet, homo_energy_bet, ave_pdos_convolved_all_bet = aux_pdos(c1, atoms, pdos_files_bet, shift, homo_occ, orbitals_cols, sigma, npoints, ave_pdos_convolved_all_bet, labels_bet) + + + ave_pdos_convolved_all_alp = np.array(ave_pdos_convolved_all_alp) + ave_pdos_convolved_total_alp = np.sum(ave_pdos_convolved_all_alp, axis=0) + + ave_pdos_convolved_all_bet = np.array(ave_pdos_convolved_all_bet) + ave_pdos_convolved_total_bet = np.sum(ave_pdos_convolved_all_bet, axis=0) - ave_pdos_convolved_all = np.array(ave_pdos_convolved_all) - ave_pdos_convolved_total = np.sum(ave_pdos_convolved_all, axis=0) + return ave_energy_grid_alp, homo_energy_alp, ave_pdos_convolved_all_alp, labels_alp, ave_pdos_convolved_total_alp, \ + ave_energy_grid_bet, homo_energy_bet, ave_pdos_convolved_all_bet, labels_bet, ave_pdos_convolved_total_bet - return ave_energy_grid, homo_energy, ave_pdos_convolved_all, labels, ave_pdos_convolved_total def exc_analysis(params): From bffa8a8d63b76b3b6fac23540ebdd401989bf88c Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Wed, 21 Feb 2024 09:21:54 -0500 Subject: [PATCH 45/47] Fix the reordering of phase gradients --- src/dyn/dyn_decoherence_methods.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index e24cca69c..43161856d 100755 --- a/src/dyn/dyn_decoherence_methods.cpp +++ b/src/dyn/dyn_decoherence_methods.cpp @@ -1540,7 +1540,7 @@ void rotate_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_par phi.set(i,i, complex(dyn_var.nab_phase[traj]->get(i, idof), 0.0) ); } } - phi = T_new.H() * phi * T_new; + phi = T_new * phi * T_new.H(); for(int i=0; i Date: Thu, 14 Mar 2024 11:18:45 -0400 Subject: [PATCH 46/47] Elaborate on the comments about the td width --- src/libra_py/dynamics/tsh/compute.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index 14eee51fd..027ffca29 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -306,11 +306,16 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): * **dyn_params["wp_width"]** ( MATRIX(ndof, 1) ): A width of a Gaussian function as an approximation to adiabatic wave packets. - This value is used for the initial conditions when the td Gaussian approximation is used, that is, use_td_width == 1 + According to the choice of the Gaussian width approximation, this parameter has different meanings + - A constant width in the fixed-width approximation, that is, `use_td_width == 0` + - The initial width in the free-particle Gaussian wave packet approximation, that is, `use_td_width == 1` + - The interaction width in the Schwarz scheme, that is, `use_td_width == 2` + - No influence on the dynamics since the width will be determined by internal variables in the Subotnik scheme, that is, `use_td_width == 3` + Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` - * **dyn_params["wp_v"]** ( MATRIX(ndof,1) ): The velocity of Gaussian wave packet following potential-free td Gaussian + * **dyn_params["wp_v"]** ( MATRIX(ndof,1) ): The velocity of Gaussian wave packet in the free-particle Gaussian approximation, that is, `use_td_width == 1` Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` @@ -342,11 +347,13 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): - 3: keep auxiliary momenta of adiabatic states except for the active state - * **dyn_params["use_td_width"]** (int): Whether to use the td Gaussian width for the nuclear wave packet approximation [ default : 0 ] - This option can be considered when it comes to unbounded systems. - This approximation is based on a nuclear wave packet on a free surface: - \sigma_x(t)=\sqrt[\sigma_x(0)^2 + (wp_v * t)^2] + * **dyn_params["use_td_width"]** (int): Options for the td Gaussian width approximations [ default : 0 ] Only used with independent-trajectory XF methods, that is, `decoherence_algo == 5 or 6` + + - 0: no td width; use the fixed-width Gaussian approximation + - 1: the td Gaussian width from a free particle Gaussian wave packet, \sigma(t)=\sqrt[\sigma(0)^2 + (wp_v * t)^2] + - 2: the Schwarz scheme where the width depends on the instantaneous de Broglie wavelength, \sigma(t)^(-2) = [\sigma(0)^2 * P/ (4 * PI) ]^2 + - 3: the Subotnik schene where the width is given pairwise depending on the auxiliary variables, \sigma_ij(t)^2 = |R_i - R_j| / |P_i - P_j| ///=============================================================================== From 35cc79cc3a996b8be4172bdafb5f3c745b3ec951 Mon Sep 17 00:00:00 2001 From: DaehoHan Date: Thu, 14 Mar 2024 11:31:37 -0400 Subject: [PATCH 47/47] Edit the explanation of the Subotnik width --- src/libra_py/dynamics/tsh/compute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index 027ffca29..32ffb2b51 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -353,7 +353,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd): - 0: no td width; use the fixed-width Gaussian approximation - 1: the td Gaussian width from a free particle Gaussian wave packet, \sigma(t)=\sqrt[\sigma(0)^2 + (wp_v * t)^2] - 2: the Schwarz scheme where the width depends on the instantaneous de Broglie wavelength, \sigma(t)^(-2) = [\sigma(0)^2 * P/ (4 * PI) ]^2 - - 3: the Subotnik schene where the width is given pairwise depending on the auxiliary variables, \sigma_ij(t)^2 = |R_i - R_j| / |P_i - P_j| + - 3: the Subotnik scheme where the width is given as a sum of pairwise widths depending on the auxiliary variables, \sigma_ij(t)^2 = |R_i - R_j| / |P_i - P_j| ///===============================================================================