diff --git a/src/dyn/Dynamics.cpp b/src/dyn/Dynamics.cpp index 6267284e2..5df9a3b20 100755 --- a/src/dyn/Dynamics.cpp +++ b/src/dyn/Dynamics.cpp @@ -570,6 +570,128 @@ 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)); + } + } + } + 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){ + 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 + } +} + +/* +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 +713,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 +735,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 +753,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 +809,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 +863,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 @@ -782,10 +915,58 @@ 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 == 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_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; + + // Propagate with the new Ham + 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 @@ -972,7 +1153,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 old_states( dyn_var.act_states); // In the interval [t, t + dt], we may have experienced the basis reordering, so we need to // change the active adiabatic state if(prms.tsh_method != 3 && prms.tsh_method != 4 ){ // Don't update states based on amplitudes, in the LZ method - dyn_var.update_active_states(); + dyn_var.update_active_states(1, 0); // 1 - forward; 0 - only active state } // For now, this function also accounts for the kinetic energy adjustments to reflect the adiabatic evolution if(prms.thermally_corrected_nbra==1){ apply_thermal_correction(dyn_var, ham, ham_aux, old_states, prms, rnd); } - // Recompute forces in respose to the updated amplitudes/density matrix/state indices + update_forces(prms, dyn_var, ham); @@ -1192,14 +1380,6 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, update_Hamiltonian_variables(prms, dyn_var, ham, ham_aux, py_funct, params, 1); - //============== Electronic propagation =================== - // Evolve electronic DOFs for all trajectories -// TEMPORARY COMMENT - //dyn_var.update_amplitudes(prms, ham); // Don't do this - then we are fine with the diabatic picture - -// dyn_var.update_density_matrix(prms, ham, 1); // This one is okay - - //============== Begin the TSH part =================== //================= Update decoherence rates & times ================ @@ -1270,44 +1450,20 @@ 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); } } // DISH, rev2023 if(prms.decoherence_algo==7){ dish_rev2023(dyn_var, ham, decoherence_rates, prms, rnd); } + + // Update amplitudes and density matrices in response to decoherence corrections dyn_var.update_amplitudes(prms, ham); dyn_var.update_density_matrix(prms, ham, 1); @@ -1360,8 +1516,17 @@ 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==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 + 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 @@ -1372,6 +1537,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); @@ -1387,7 +1553,9 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params, else{ cout<<"tsh_method == "<n_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; project_out_aux = x.project_out_aux; + tp_algo = x.tp_algo; + use_td_width = x.use_td_width; ///================= Entanglement of trajectories ================================ entanglement_opt = x.entanglement_opt; @@ -184,6 +196,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); @@ -206,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; @@ -215,7 +230,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{ @@ -343,10 +359,26 @@ 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_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]); } 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]); } @@ -378,6 +410,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..90152b9a1 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 @@ -389,6 +390,7 @@ class dyn_control_params{ - 5: SHXF of Min - 6: MQCXF - 7: DISH, rev2023 + - 8: diabatic IDA, experimental */ double decoherence_algo; @@ -464,8 +466,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; @@ -496,10 +503,18 @@ class dyn_control_params{ /** - The width of frozen 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; + + + /** + MATRIX(ndof, 1) of wave packet velocities for the decoherence from SHXF & MQCXF + This value is applied when use_td_width = 1 + [ default : NULL ] + */ + MATRIX* wp_v; /** @@ -507,7 +522,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 @@ -515,12 +538,36 @@ 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 + + Options: + - 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 + */ + 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; + + /** 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. @@ -695,6 +742,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 diff --git a/src/dyn/dyn_decoherence.h b/src/dyn/dyn_decoherence.h index b3d4f7b00..164332ca6 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); @@ -72,13 +75,15 @@ CMATRIX mfsd(MATRIX& p, CMATRIX& Coeff, MATRIX& invM, double dt, vector& // 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 // XF propagation +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, int do_rotation); +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); ///================ In dyn_decoherence_time.cpp =================================== diff --git a/src/dyn/dyn_decoherence_methods.cpp b/src/dyn/dyn_decoherence_methods.cpp index f24dd101f..43161856d 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 = = "<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; } } @@ -904,7 +1038,9 @@ void xf_destroy_AT(dyn_variables& dyn_var, nHamiltonian& ham, double threshold){ double alpha = compute_kinetic_energy(p_real, *dyn_var.iM) + Epot_old - Epot; if(alpha > 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)); @@ -939,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; @@ -988,6 +1127,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]; @@ -996,6 +1136,8 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn for(int i=0; iget(idof, traj));} @@ -1018,11 +1160,13 @@ 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]; int a = dyn_var.act_states[traj]; - + MATRIX p_real(ndof, 1); MATRIX p_aux_temp(ndof, 1); @@ -1037,9 +1181,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){ @@ -1070,19 +1217,60 @@ void shxf(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn } // Check the turning point - if (is_first[i] == 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; 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));} @@ -1155,8 +1348,12 @@ 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]; + + int a = dyn_var.act_states[traj]; MATRIX p_real(ndof, 1); MATRIX p_aux_temp(ndof, 1); @@ -1175,23 +1372,26 @@ 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){ 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; - 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 <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; - - 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; - } + if (is_first[i] == 0 and prms.tp_algo != 0){ + MATRIX Fa(ndof, 1); + for(int idof=0; idofget(idof,0) ); + } + double dp_old = 0.0; + double dp_new = 0.0; - for(int idof=0; idofset(idof, traj, dyn_var.p->get(idof, traj) * sqrt(alpha)); - } + 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 << "Total 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 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 } // traj } -void XF_correction(CMATRIX& Ham, dyn_variables& dyn_var, CMATRIX& C, double wp_width, CMATRIX& T, int traj){ +void rotate_nab_phase(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms){ + int ntraj = dyn_var.ntraj; + 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]; + + CMATRIX T_new(*dyn_var.proj_adi[traj]); + if(prms.assume_always_consistent){ T_new.identity(); } + + // 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) ); + } + } + phi = T_new * phi * T_new.H(); + + for(int i=0; iset(i, idof, phi.get(i,i).real() ); + } + } + }//idof + + if (prms.decoherence_algo==5){ + for(int i=0; iset(i, idof, dyn_var.nab_phase[traj]->get(i, idof) ); + } + } + } + } + + } // 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 = T * C * C.H() * T.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(wp_width, 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)) ); + if(is_mixed[i]==1){ + 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))); + } + } } - F = T * F * T.H(); - 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){ @@ -1332,9 +1631,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) ); + if(is_mixed[i]==1){ + 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 *dyn_var.f += *dyn_var.f_xf; } -void propagate_half_xf(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms, int do_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; @@ -1411,6 +1704,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=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(); - - // Don't apply T, since hamiltonian elements were already transformed through the transform_all method - CMATRIX T_new(nadi, nadi); T_new.load_identity(); - //P = ham->get_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 ); - // - //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) )] - - XF_correction(Hxf_old, dyn_var, C, prms.wp_width, T, itraj); - XF_correction(Hxf, dyn_var, C, prms.wp_width, T_new, itraj); - - Hxf = T_new.H() * Hxf * T_new; - Hxf += Hxf_old; - - 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; - } + // Generate the half-time exponential operator + Hxf = *dyn_var.ham_xf[itraj]; + D = libspecialfunctions::exp_(Hxf, complex(0.0, -0.5*dt) ); - *dyn_var.proj_adi[itraj] = T_new; + C = D * C; // Insert the propagated result back for(int st=0; stget_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; + + double diff = 0.0; + + 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(); + + 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 + + + }// 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 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); diff --git a/src/dyn/dyn_variables.cpp b/src/dyn/dyn_variables.cpp index 4a91964dc..38bfffbae 100755 --- a/src/dyn/dyn_variables.cpp +++ b/src/dyn/dyn_variables.cpp @@ -181,9 +181,13 @@ void dyn_variables::allocate_shxf(){ for(int itraj=0; itraj()); is_first.push_back(vector()); + is_fixed.push_back(vector()); + is_keep.push_back(vector()); for(int i=0; i(ntraj); p_aux_old = vector(ntraj); nab_phase = vector(ntraj); + nab_phase_old = vector(ntraj); + ham_xf = vector(ntraj); for(int itraj=0; itraj()); is_first.push_back(vector()); + is_fixed.push_back(vector()); + is_keep.push_back(vector()); for(int i=0; i(ntraj); p_aux_old = vector(ntraj); nab_phase = vector(ntraj); + //nab_phase_old = vector(ntraj); + ham_xf = vector(ntraj); for(int itraj=0; itraj > 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 @@ -372,6 +388,30 @@ class dyn_variables{ vector */ vector nab_phase; + + /** + Spatial derivative of the phase of coefficients of state-wise auxiliary trajectories + + Options: + vector + */ + vector nab_phase_old; + + /** + XF Hamiltonian + + Options: + vector + */ + vector ham_xf; + + /** + Wave packet widths based on the Gaussian approximation + + Options: + MATRIX(ndof, ntraj) + */ + MATRIX* wp_width; /** Quantum momenta defined as (-1) * \nabla_nuc |\chi| / |\chi| @@ -425,6 +465,13 @@ class dyn_variables{ vector tcnbra_ekin; + ///================= Misc =================== + /** + The current MD time step + */ + int timestep; + + ///====================== In dyn_variables.cpp ===================== @@ -460,6 +507,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; } @@ -467,7 +515,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;} + } + } @@ -494,6 +549,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); 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; diff --git a/src/dyn/libdyn.cpp b/src/dyn/libdyn.cpp index 1160f2572..af74caf00 100755 --- a/src/dyn/libdyn.cpp +++ b/src/dyn/libdyn.cpp @@ -103,8 +103,10 @@ 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) ///================= Entanglement of trajectories ================================ .def_readwrite("entanglement_opt", &dyn_control_params::entanglement_opt) @@ -122,6 +124,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 +194,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__) @@ -241,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) @@ -267,7 +274,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) @@ -330,6 +338,10 @@ void export_dyn_decoherence_objects(){ int instantaneous_decoherence_variant, int collapse_option) = &instantaneous_decoherence; def("instantaneous_decoherence", expt_instantaneous_decoherence_v1); + void (*expt_instantaneous_decoherence_dia_v1)(CMATRIX& Coeff, nHamiltonian& ham, + vector& 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) @@ -354,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; @@ -705,6 +721,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; 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/__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", diff --git a/src/libra_py/data_conv.py b/src/libra_py/data_conv.py index a2ddc8901..34bc780cd 100755 --- a/src/libra_py/data_conv.py +++ b/src/libra_py/data_conv.py @@ -578,3 +578,4 @@ def vasp_to_xyz(filename): f.close() + 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/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: diff --git a/src/libra_py/dynamics/tsh/compute.py b/src/libra_py/dynamics/tsh/compute.py index 750ddf9d3..32ffb2b51 100755 --- a/src/libra_py/dynamics/tsh/compute.py +++ b/src/libra_py/dynamics/tsh/compute.py @@ -305,7 +305,17 @@ 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. + 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 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` @@ -313,15 +323,39 @@ 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 Ha ] + 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` - 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): 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 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| + + ///=============================================================================== ///================= Entanglement of trajectories ================================ ///=============================================================================== @@ -514,7 +548,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 @@ -534,8 +568,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 + "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 } ) #================= Entanglement of trajectories ================================ 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") 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)) 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..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): @@ -2020,7 +2105,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): """ 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