Skip to content

Implement ITMQC-XF methods in JCP, 156, 174109 (2022) #228

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jkha-unist wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Implement ITMQC-XF methods in JCP, 156, 174109 (2022) #228

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
29dbae8
Copy raw shxf.py used for the paper (based on old master: 3e273fa)
jkha-unist Jan 17, 2023
3b7c1dc
Remove XF1D
jkha-unist Jan 18, 2023
f980680
Sync and trim variables in shxf.py with current master
jkha-unist Jan 25, 2023
cd7318a
Remove istep in hopping-related routines and add resetting l_fix vari…
jkha-unist Jan 25, 2023
f306ac5
Change the order of if-branch for aux_propagator
jkha-unist Jan 25, 2023
3bb4a5d
Test comment recognition
jkha-unist Jan 25, 2023
dc87661
Add 'fix' functionality in aux_propagator
jkha-unist Jan 25, 2023
1539828
Add k_lk in el_prop/*xf*
jkha-unist Jan 27, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
162 changes: 28 additions & 134 deletions src/mqc/el_prop/derivs_xf.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,174 +4,68 @@
#include <math.h>

// Routine to calculate cdot contribution originated from XF term
static void xf_cdot(int nat, int ndim, int nst, int *l_coh, double *mass, double *sigma,
double **pos, double **qmom, double ***aux_pos, double ***phase, double complex *c, double complex *xfcdot){
// CAUTION!!! different from the k_lk in CT!!!
static void xf_cdot(int nst, double **k_lk, double complex *c, double complex *xfcdot){

double **dec = malloc(nst * sizeof(double*));
double *rho = malloc(nst * sizeof(double));

int ist, jst, iat, isp;

// Initialize variables related to decoherence
for(iat = 0; iat < nat; iat++){
for(isp = 0; isp < ndim; isp++){
qmom[iat][isp] = 0.0;
}
}

for(ist = 0; ist < nst; ist++){
dec[ist] = malloc(nst * sizeof(double));
for(jst = 0; jst < nst; jst++){
dec[ist][jst] = 0.0;
}
}
int lst, mst;

// Calculate densities from current coefficients
for(ist = 0; ist < nst; ist++){
rho[ist] = creal(conj(c[ist]) * c[ist]);
}

// Get quantum momentum from auxiliary positions and sigma values
for(ist = 0; ist < nst; ist++){

if(l_coh[ist] == 1){
for(iat = 0; iat < nat; iat++){
for(isp = 0; isp < ndim; isp++){
qmom[iat][isp] += 0.5 * rho[ist] * (pos[iat][isp] - aux_pos[ist][iat][isp])
/ pow(sigma[iat], 2.0) / mass[iat];
}
}
}

for(lst = 0; lst < nst; lst++){
rho[lst] = creal(conj(c[lst]) * c[lst]);
}

// Get decoherence term from quantum momentum and phase
for(ist = 0; ist < nst; ist++){
for(jst = ist + 1; jst < nst; jst++){

if(l_coh[ist] == 1 && l_coh[jst] == 1){
for(iat = 0; iat < nat; iat++){
for(isp = 0; isp < ndim; isp++){
dec[ist][jst] += qmom[iat][isp] * (phase[ist][iat][isp] - phase[jst][iat][isp]);
}
}
}
dec[jst][ist] = - 1.0 * dec[ist][jst];

// Get cdot contribution from decoherence term, state-pair expression
for(lst = 0; lst < nst; lst++){
xfcdot[lst] = 0.0 + 0.0 * I;
for(mst = 0; mst < nst; mst++){
xfcdot[lst] -= rho[mst] * k_lk[mst][lst] * c[lst];
}
}

// Get cdot contribution from decoherence term
for(ist = 0; ist < nst; ist++){
xfcdot[ist] = 0.0 + 0.0 * I;
for(jst = 0; jst < nst; jst++){
xfcdot[ist] -= rho[jst] * dec[jst][ist] * c[ist];
}
}

// Deallocate temporary arrays
for(ist = 0; ist < nst; ist++){
free(dec[ist]);
}

free(dec);
free(rho);

}

// Routine to print xf debug info
static void xf_print_coef(int nst, double complex *coef, double complex *xfcdot, double *dotpopdec){
int ist;
int lst;

for(ist = 0; ist < nst; ist++){
dotpopdec[ist] = 2.0 * creal(xfcdot[ist] * conj(coef[ist]));
for(lst = 0; lst < nst; lst++){
dotpopdec[lst] = 2.0 * creal(xfcdot[lst] * conj(coef[lst]));
}
}

// CAUTION!!! different from the k_lk in CT!!!
// Routine to calculate rhodot contribution originated from XF term
static void xf_rhodot(int nat, int ndim, int nst, int *l_coh, double *mass, double *sigma,
double **pos, double **qmom, double ***aux_pos, double ***phase, double complex **rho, double complex **xfrhodot){

double **dec = malloc(nst * sizeof(double*));
static void xf_rhodot(int nst, double **k_lk, double complex **rho, double complex **xfrhodot){

int ist, jst, kst, iat, isp;

// Initialize variables related to decoherence
for(iat = 0; iat < nat; iat++){
for(isp = 0; isp < ndim; isp++){
qmom[iat][isp] = 0.0;
}
}

for(ist = 0; ist < nst; ist++){
dec[ist] = malloc(nst * sizeof(double));
for(jst = 0; jst < nst; jst++){
dec[ist][jst] = 0.0;
}
}

// Get quantum momentum from auxiliary positions and sigma values
for(ist = 0; ist < nst; ist++){

if(l_coh[ist] == 1){
for(iat = 0; iat < nat; iat++){
for(isp = 0; isp < ndim; isp++){
qmom[iat][isp] += 0.5 * creal(rho[ist][ist]) * (pos[iat][isp] - aux_pos[ist][iat][isp])
/ pow(sigma[iat], 2.0) / mass[iat];
}
}
}

}

// Get decoherence term from quantum momentum and phase
for(ist = 0; ist < nst; ist++){
for(jst = ist + 1; jst < nst; jst++){

if(l_coh[ist] == 1 && l_coh[jst] == 1){
for(iat = 0; iat < nat; iat++){
for(isp = 0; isp < ndim; isp++){
dec[ist][jst] += qmom[iat][isp] * (phase[ist][iat][isp] - phase[jst][iat][isp]);
}
}
}
dec[jst][ist] = - 1.0 * dec[ist][jst];

}
}
int lst, kst, mst;

// Get rhodot contribution from decoherence term
for(ist = 0; ist < nst; ist++){
for(lst = 0; lst < nst; lst++){
// Diagonal components
xfrhodot[ist][ist] = 0.0 + 0.0 * I;
for(kst = 0; kst < nst; kst++){
xfrhodot[ist][ist] -= 2.0 * dec[kst][ist] * rho[ist][kst] * rho[kst][ist];
xfrhodot[lst][lst] = 0.0 + 0.0 * I;
for(mst = 0; mst < nst; mst++){
xfrhodot[lst][lst] -= 2.0 * k_lk[mst][lst] * rho[lst][mst] * rho[mst][lst];
}
// Off-diagonal components
for(jst = ist + 1; jst < nst; jst++){
xfrhodot[ist][jst] = 0.0 + 0.0 * I;
for(kst = 0; kst < nst; kst++){
xfrhodot[ist][jst] -= (dec[kst][ist] + dec[kst][jst]) * rho[ist][kst] * rho[kst][jst];
for(kst = lst + 1; kst < nst; kst++){
xfrhodot[lst][kst] = 0.0 + 0.0 * I;
for(mst = 0; mst < nst; mst++){
xfrhodot[lst][kst] -= (k_lk[mst][lst] + k_lk[mst][kst]) * rho[lst][mst] * rho[mst][kst];
}
xfrhodot[jst][ist] = conj(xfrhodot[ist][jst]);
xfrhodot[kst][lst] = conj(xfrhodot[lst][kst]);
}
}

// Deallocate temporary arrays
for(ist = 0; ist < nst; ist++){
free(dec[ist]);
}

free(dec);

}

// Routine to print xf debug info
static void xf_print_rho(int nst, double complex **xfrhodot, double *dotpopdec){
int ist;
int lst;

for(ist = 0; ist < nst; ist++){
dotpopdec[ist] = creal(xfrhodot[ist][ist]);
for(lst = 0; lst < nst; lst++){
dotpopdec[lst] = creal(xfrhodot[lst][lst]);
}
}

86 changes: 11 additions & 75 deletions src/mqc/el_prop/el_propagator_xf.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,69 +5,43 @@ import numpy as np
cimport numpy as np

cdef extern from "rk4_xf.c":
void rk4(int nat, int ndim, int nst, int nesteps, double dt, char *elec_object, \
bint *l_coh, double *mass, double *energy, double *energy_old, double *sigma, \
double **nacme, double **nacme_old, double **pos, double **qmom, double ***aux_pos, \
double ***phase, double complex *coef, double complex **rho, int verbosity, \
void rk4(int nst, int nesteps, double dt, char *elec_object, double *energy, double *energy_old,
double **nacme, double **nacme_old, double **k_lk, double complex *coef, double complex **rho, int verbosity, \
double *dotpopdec)

def el_run(md):
cdef:
char *elec_object_c
bint *l_coh
double *mass
double *energy
double *energy_old
double *sigma
double **nacme
double **nacme_old
double **pos
double **qmom
double ***aux_pos
double ***phase
double **k_lk
double complex *coef
double complex **rho
double *dotpopdec

bytes py_bytes
int ist, jst, nst, nesteps, iat, aux_nat, aux_ndim, verbosity
int ist, jst, nst, nesteps, verbosity
double dt

# Assign size variables
nst = md.mol.nst
nesteps, dt = md.nesteps, md.dt
aux_nat, aux_ndim = md.aux.nat, md.aux.ndim

# Allocate variables
l_coh = <bint*> PyMem_Malloc(nst * sizeof(bint))
mass = <double*> PyMem_Malloc(aux_nat * sizeof(double))
energy = <double*> PyMem_Malloc(nst * sizeof(double))
energy_old = <double*> PyMem_Malloc(nst * sizeof(double))
sigma = <double*> PyMem_Malloc(aux_nat * sizeof(double))

nacme = <double**> PyMem_Malloc(nst * sizeof(double*))
nacme_old = <double**> PyMem_Malloc(nst * sizeof(double*))
pos = <double**> PyMem_Malloc(aux_nat * sizeof(double*))
qmom = <double**> PyMem_Malloc(aux_nat * sizeof(double*))

aux_pos = <double***> PyMem_Malloc(nst * sizeof(double**))
phase = <double***> PyMem_Malloc(nst * sizeof(double**))
k_lk = <double**> PyMem_Malloc(nst * sizeof(double*))

for ist in range(nst):
nacme[ist] = <double*> PyMem_Malloc(nst * sizeof(double))
nacme_old[ist] = <double*> PyMem_Malloc(nst * sizeof(double))
k_lk[ist] = <double*> PyMem_Malloc(nst * sizeof(double))

for iat in range(aux_nat):
pos[iat] = <double*> PyMem_Malloc(aux_ndim * sizeof(double))
qmom[iat] = <double*> PyMem_Malloc(aux_ndim * sizeof(double))

for ist in range(nst):
aux_pos[ist] = <double**> PyMem_Malloc(aux_nat * sizeof(double*))
phase[ist] = <double**> PyMem_Malloc(aux_nat * sizeof(double*))
for iat in range(aux_nat):
aux_pos[ist][iat] = <double*> PyMem_Malloc(aux_ndim * sizeof(double))
phase[ist][iat] = <double*> PyMem_Malloc(aux_ndim * sizeof(double))

# Debug related
verbosity = md.verbosity
if (verbosity >= 1):
Expand All @@ -78,25 +52,11 @@ def el_run(md):
energy[ist] = md.mol.states[ist].energy
energy_old[ist] = md.mol.states[ist].energy_old

for iat in range(aux_nat):
mass[iat] = md.aux.mass[iat]
sigma[iat] = md.sigma[iat]

for ist in range(nst):
for jst in range(nst):
nacme[ist][jst] = md.mol.nacme[ist, jst]
nacme_old[ist][jst] = md.mol.nacme_old[ist, jst]

for iat in range(aux_nat):
for isp in range(aux_ndim):
pos[iat][isp] = md.pos_0[iat, isp]

for ist in range(nst):
l_coh[ist] = md.l_coh[ist]
for iat in range(aux_nat):
for isp in range(aux_ndim):
aux_pos[ist][iat][isp] = md.aux.pos[ist, iat, isp]
phase[ist][iat][isp] = md.phase[ist, iat, isp]
k_lk[ist][jst] = md.k_lk[ist, jst]

# Assign coef or rho with respect to propagation scheme
if (md.elec_object == "coefficient"):
Expand All @@ -121,8 +81,8 @@ def el_run(md):

# Propagate electrons depending on the propagator
if (md.propagator == "rk4"):
rk4(aux_nat, aux_ndim, nst, nesteps, dt, elec_object_c, l_coh, mass, energy, \
energy_old, sigma, nacme, nacme_old, pos, qmom, aux_pos, phase, coef, rho, verbosity, dotpopdec)
rk4(nst, nesteps, dt, elec_object_c, energy, energy_old, nacme, nacme_old, \
k_lk, coef, rho, verbosity, dotpopdec)

# Assign variables from C to python
if (md.elec_object == "coefficient"):
Expand Down Expand Up @@ -155,42 +115,18 @@ def el_run(md):
if (jst != ist):
md.dotpopnac[ist] -= 2. * nacme[ist][jst] * md.mol.rho.real[jst, ist]

if (verbosity >= 2):
for iat in range(aux_nat):
for isp in range(aux_ndim):
md.qmom[iat, isp] = qmom[iat][isp]

# Deallocate variables
for ist in range(nst):
for iat in range(aux_nat):
PyMem_Free(aux_pos[ist][iat])
PyMem_Free(phase[ist][iat])

for ist in range(nst):
PyMem_Free(nacme[ist])
PyMem_Free(nacme_old[ist])
PyMem_Free(k_lk[ist])

for iat in range(aux_nat):
PyMem_Free(pos[iat])
PyMem_Free(qmom[iat])

for ist in range(nst):
PyMem_Free(aux_pos[ist])
PyMem_Free(phase[ist])

PyMem_Free(l_coh)
PyMem_Free(mass)
PyMem_Free(energy)
PyMem_Free(energy_old)
PyMem_Free(sigma)

PyMem_Free(nacme)
PyMem_Free(nacme_old)
PyMem_Free(pos)
PyMem_Free(qmom)

PyMem_Free(aux_pos)
PyMem_Free(phase)
PyMem_Free(k_lk)

if (verbosity >= 1):
PyMem_Free(dotpopdec)
Expand Down
Loading