Skip to content

Commit 1826c4c

Browse files
alexvakimovalexvakimov
authored
Merge pull request #252 from suarez-va/devel
Restructured code to minimize number of function calls, added in KCRPMD forces. Created new RPMD file for general RPMD calls
2 parents 69b2cc1 + b9ee0be commit 1826c4c

9 files changed

Lines changed: 431 additions & 142 deletions

File tree

‎src/dyn/dyn_control_params.cpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,9 @@ dyn_control_params::dyn_control_params(){
8484
use_qtsh = 0;
8585
qtsh_force_option = 1;
8686

87+
///================= KCRPMD specific ====================
88+
use_kcrpmd = 0;
89+
8790
///================= Decoherence options =========================================
8891
decoherence_algo = -1;
8992
sdm_norm_tolerance = 0.0;

‎src/dyn/dyn_control_params.h‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -58,6 +58,7 @@ class dyn_control_params{
5858
- 1: adiabatic representation, wfc [ default ]
5959
- 2: diabatic representation, density matrix (e.g. Liouville's picture)
6060
- 3: adiabatic representation, density matrix (e.g. Liouville's picture)
61+
- 4: MMST adiabatic representation, wfc
6162
*/
6263
int rep_tdse;
6364

@@ -132,6 +133,7 @@ class dyn_control_params{
132133
- 1: state-specific as in the TSH or adiabatic (including adiabatic excited states) [ default ]
133134
- 2: Ehrenfest and MMST (SQC)
134135
- 3: QTSH force
136+
- 4: KC-RPMD force
135137
*/
136138
int force_method;
137139

@@ -493,6 +495,18 @@ class dyn_control_params{
493495
*/
494496
int qtsh_force_option;
495497

498+
499+
//=========== KC-RPMD options ==========
500+
/**
501+
Whether to use KC-RPMD
502+
503+
Options:
504+
505+
- 0: don't apply [ default ]
506+
- 1: use it
507+
*/
508+
int use_kcrpmd;
509+
496510

497511
///===============================================================================
498512
///================= Decoherence options =========================================
@@ -756,6 +770,7 @@ class dyn_control_params{
756770
- 1: ETHD
757771
- 2: ETHD3 (experimental)
758772
- 22: another flavor of ETHD3 (experimental)
773+
- 3: RPMD
759774
*/
760775
int entanglement_opt;
761776

‎src/dyn/dyn_variables.cpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -310,6 +310,8 @@ void dyn_variables::allocate_kcrpmd(){
310310
if(kcrpmd_vars_status==0){
311311

312312
auxiliary_y = vector<double>(1, -1.0);
313+
auxiliary_vy = vector<double>(1, 0.0);
314+
auxiliary_fy = vector<double>(1, 0.0);
313315

314316
kcrpmd_vars_status = 1;
315317
}
@@ -482,6 +484,8 @@ dyn_variables::dyn_variables(const dyn_variables& x){
482484

483485
// Copy content
484486
auxiliary_y = x.auxiliary_y;
487+
auxiliary_vy = x.auxiliary_vy;
488+
auxiliary_fy = x.auxiliary_fy;
485489

486490
}// if KCRPMD vars
487491

@@ -644,6 +648,8 @@ dyn_variables::~dyn_variables(){
644648

645649
if(kcrpmd_vars_status==1){
646650
auxiliary_y.clear();
651+
auxiliary_vy.clear();
652+
auxiliary_fy.clear();
647653

648654
kcrpmd_vars_status = 0;
649655
}

‎src/dyn/dyn_variables.h‎

Lines changed: 11 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -529,7 +529,7 @@ class dyn_variables{
529529
MATRIX* qtsh_f_nc;
530530

531531

532-
///========= For kinetic-constrained RPMD ======================
532+
///========= For KC-RPMD ======================
533533
/**
534534
Status of the KC-RPMD vars
535535
@@ -539,10 +539,19 @@ class dyn_variables{
539539
int kcrpmd_vars_status;
540540

541541
/**
542-
The classical auxiliary electronic variable(s)
542+
The classical auxiliary electronic variable
543543
*/
544544
vector<double> auxiliary_y;
545545

546+
/**
547+
The classical auxiliary electronic variable velocity
548+
*/
549+
vector<double> auxiliary_vy;
550+
551+
/**
552+
The classical auxiliary electronic variable force
553+
*/
554+
vector<double> auxiliary_fy;
546555

547556

548557
///================= Misc ===================

‎src/libra_py/dynamics/tsh/compute.py‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,13 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd):
126126
- 1: state-specific as in the TSH or adiabatic (including adiabatic excited states) [ default ]
127127
- 2: Ehrenfest
128128
- 3: QTSH force
129+
- 4: KC-RPMD force
130+
131+
132+
* **dyn_params["sqc_gamma"]** ( double ): Gamma parameter in the Meyer-Miller Hamiltonian of form:
133+
H = sum_{i,j} H_{ij} [ 1/(2 hbar) * (q_i - i * p_j) * (q_j + i * p_i) - gamma * delta_ij ]
134+
135+
- 0.0: corresponds to Ehrenfest force [ default ]
129136
130137
131138
* **dyn_params["enforce_state_following"]** ( int ): Wheather we want to enforce nuclear dynamics
@@ -203,7 +210,8 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd):
203210
- -1: use LD approach, it includes phase correction too [ default ]
204211
- 0: no state tracking
205212
- 1: method of Kosuke Sato (may fail by getting trapped into an infinite loop)
206-
- 2: Munkres-Kuhn (Hungarian) algorithm
213+
- 2: Munkres-Kuhn (Hungarian) algorithm
214+
- 21: ChatGPT-generated Munkres-Kuhn (Hungarian) algorithm
207215
- 3: experimental stochastic algorithm, the original version with elimination (known problems)
208216
- 32: experimental stochastic algorithms with all permutations (too expensive)
209217
- 33: the improved stochastic algorithm with good scaling and performance, on par with the mincost
@@ -255,6 +263,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd):
255263
- 6: MASH
256264
- 7: FSSH2
257265
- 8: FSSH3
266+
- 9: GFSH (original)
258267
259268
260269
* **dyn_params["hop_acceptance_algo"]** ( int ): Options to control the acceptance of the proposed hops:
@@ -368,6 +377,15 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd):
368377
- 0: Considering only the first-order force, i.e., off-diagonal Ehrenfest force
369378
- 1: The whole force including the second-order term is used [default]
370379
380+
381+
//=========== KC-RPMD options ==========
382+
383+
* **dyn_params["use_kcrpmd"]** ( int ): Whether to use KC-RPMD
384+
385+
- 0: don't apply [ default ]
386+
- 1: use it
387+
388+
371389
///===============================================================================
372390
///================= Decoherence options =========================================
373391
///===============================================================================
@@ -549,6 +567,7 @@ def run_dynamics(dyn_var, _dyn_params, ham, compute_model, _model_params, rnd):
549567
- 1: ETHD
550568
- 2: ETHD3 (experimental)
551569
- 22: another flavor of ETHD3 (experimental)
570+
- 3: RPMD
552571
553572
554573
* **dyn_params["ETHD3_alpha"]** ( double ): Gaussian exponents that dresses up the trajectories in the ETHD3 method

‎src/nhamiltonian/libnhamiltonian.cpp‎

Lines changed: 5 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -317,12 +317,9 @@ void export_nhamiltonian_objects(){
317317
// = &nHamiltonian::QTSH_forces_dia;
318318

319319

320-
double (nHamiltonian::*expt_internal_ring_polymer_potential_v1)(const MATRIX& q, const MATRIX& invM, double beta) = &nHamiltonian::internal_ring_polymer_potential;
321-
vector<MATRIX> (nHamiltonian::*expt_generate_m_matrices_v1)(double beta) = &nHamiltonian::generate_m_matrices;
322-
double (nHamiltonian::*expt_kinked_pair_electronic_potential_v1)(double beta) = &nHamiltonian::kinked_pair_electronic_potential;
323-
double (nHamiltonian::*expt_kinetic_constraint_potential_v1)(double beta, double eta, double b) = &nHamiltonian::kinetic_constraint_potential;
324-
double (nHamiltonian::*expt_heavy_side_potential_v1)(vector<double>& auxiliary_y, int theta, double beta, double b) = &nHamiltonian::heavy_side_potential;
325-
double (nHamiltonian::*expt_kcrpmd_effective_potential_v1)(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b) = &nHamiltonian::kcrpmd_effective_potential;
320+
double (nHamiltonian::*expt_kcrpmd_effective_potential_v1)(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b, double c, double d) = &nHamiltonian::kcrpmd_effective_potential;
321+
MATRIX (nHamiltonian::*expt_kcrpmd_effective_force_v1)(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b, double c, double d) = &nHamiltonian::kcrpmd_effective_force;
322+
vector<double> (nHamiltonian::*expt_kcrpmd_effective_auxiliary_force_v1)(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b, double c, double d) = &nHamiltonian::kcrpmd_effective_auxiliary_force;
326323

327324

328325
void (nHamiltonian::*expt_add_ethd_dia_v1)(const MATRIX& q, const MATRIX& invM, int der_lvl) = &nHamiltonian::add_ethd_dia;
@@ -592,12 +589,9 @@ void export_nhamiltonian_objects(){
592589
// .def("QTSH_forces_dia", expt_QTSH_forces_dia_v2)
593590

594591

595-
.def("internal_ring_polymer_potential", expt_internal_ring_polymer_potential_v1)
596-
.def("generate_m_matrices", expt_generate_m_matrices_v1)
597-
.def("kinked_pair_electronic_potential", expt_kinked_pair_electronic_potential_v1)
598-
.def("kinetic_constraint_potential", expt_kinetic_constraint_potential_v1)
599-
.def("heavy_side_potential", expt_heavy_side_potential_v1)
600592
.def("kcrpmd_effective_potential", expt_kcrpmd_effective_potential_v1)
593+
.def("kcrpmd_effective_force", expt_kcrpmd_effective_force_v1)
594+
.def("kcrpmd_effective_auxiliary_force", expt_kcrpmd_effective_auxiliary_force_v1)
601595

602596
;
603597

‎src/nhamiltonian/nHamiltonian.h‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -598,14 +598,10 @@ class nHamiltonian{
598598

599599

600600
///< In nHamiltonian_compute_KCRPMD.cpp
601-
double internal_ring_polymer_potential(const MATRIX& q, const MATRIX& invM, double beta);
602601
vector<MATRIX> generate_m_matrices(double beta);
603-
double donor_electronic_potential(double beta);
604-
double acceptor_electronic_potential(double beta);
605-
double kinked_pair_electronic_potential(double beta);
606-
double kinetic_constraint_potential(double beta, double eta, double a);
607-
double heavy_side_potential(vector<double>& auxiliary_y, int theta, double beta, double b);
608-
double kcrpmd_effective_potential(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b);
602+
double kcrpmd_effective_potential(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b, double c, double d);
603+
MATRIX kcrpmd_effective_force(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b, double c, double d);
604+
vector<double> kcrpmd_effective_auxiliary_force(vector<double>& auxiliary_y, const MATRIX& q, const MATRIX& invM, double beta, double eta, double a, double b, double c, double d);
609605

610606

611607

@@ -645,6 +641,11 @@ MATRIX ETHD3_forces(const MATRIX& q, const MATRIX& p, const MATRIX& invM, double
645641
MATRIX ETHD3_friction(const MATRIX& q, const MATRIX& p, const MATRIX& invM, double alp, double bet); // dH_{ETHD3}/dP
646642

647643

644+
///< In nHamiltonian_compute_RPMD.cpp
645+
double RPMD_internal_potential(const MATRIX& q, const MATRIX& invM, double beta);
646+
MATRIX RPMD_internal_force(const MATRIX& q, const MATRIX& invM, double beta);
647+
648+
648649

649650

650651
}// namespace libnhamiltonian

0 commit comments

Comments
 (0)