Conversation
Add FSSH using kTDC method
|
@felixmiaomiao @puzhichen Please review the refactored code and APIs and let me know if they make sense. |
| import unittest | ||
| import numpy as np | ||
| import cupy as cp | ||
| import scipy.linalg |
| # mo2_reordered[:,i] *= -1 | ||
| sign_array[i] = -1 | ||
|
|
||
| sign_array = np.sign(cp.asnumpy(final_chosen_overlaps)) |
There was a problem hiding this comment.
(match_and_reorder_mos is not used, so the code may leave unchanged.)
Should we consider the reordering of MO coefficients here? In the original version, s_mo_new was re-calculated to ensure the phase (sign) of the reordered orbitals remains consistent with the reference. Since final_chosen_overlaps is derived from absolute values, using np.sign on it will always return 1, which might skip the necessary phase correction. Even if the new _wfn_overlap function no longer relies on mo2_reordered, leaving this broken sign_array logic here might cause hidden bugs if this function is reused elsewhere.
| wvo, | ||
| max_cycle=td_nac.cphf_max_cycle, | ||
| tol=td_nac.cphf_conv_tol)[0] # eq.(80) in Ref. [1] | ||
| z1 /= EJ-EI # only one spin, negative in cphf |
There was a problem hiding this comment.
This modification is critical when calculating combinations of multiple states. I'm concerned about the numerical stability of the two approaches—I should verify which one is more robust in these scenarios.
| if states[0] != 0: # excited state | ||
| sign = check_phase_modified(mol0, mo_coeff0, mo2_reordered, xi0, xi1, nocc, s) | ||
| self.sign *= np.sign(sign) | ||
| state_ovlp = _wfn_overlap(mo_coeff0, mo_coeff, xi0, xi1, s) |
There was a problem hiding this comment.
I noticed we are passing the un-reordered mo_coeff here instead of mo2_reordered. I assume this is intentional because the excitation amplitudes xi1 are in the basis of the original mo_coeff, and the Plasser algorithm (_wfn_overlap) naturally accounts for any orbital permutations and phase changes via the determinant of the raw overlap matrix. Could you confirm if my understanding is correct?
| std_dev = np.sqrt(k_b * temperature / masses_amu[:, np.newaxis]) | ||
| velocities = np.random.normal(0, std_dev, (N, 3)) | ||
| if force_temp: | ||
| total_momentum = np.sum(velocities * masses_amu[:, np.newaxis], axis=0) |
There was a problem hiding this comment.
Could we consider moving the center-of-mass momentum removal outside the if force_temp block to prevent the molecule from drifting? I understand our current code always uses force_temp=True, but this would make the function more robust if the default False is ever used in the future.
| nvir = nmo - nocc | ||
| if isinstance(td_scanner, RisBase): | ||
| xs0 = td_scanner.xy[0] | ||
| xs0 = [xs0[i-1].reshape(nocc, nvir) for i in self.states] |
There was a problem hiding this comment.
If ground state index 0 is in self.states, will this part be problematic?
|
|
||
| states_reorder = [] | ||
| for i in self.states: | ||
| state_ovlp = _wfn_overlap(mo_coeff0, mo_coeff, xs0[i-1], xs1[i-1], s) |
There was a problem hiding this comment.
Can self.states be non-contiguous? The lengths of xs0 and xs1 are len(self.states). Indexin with i-1 here may be out of bounds.
|
|
||
| # Run a FSSH simulation including only the first and second excited | ||
| # electronic states. This restricts surface hopping among these states. | ||
| fssh = FSSH_TDDFT(td, [1, 2]) |
There was a problem hiding this comment.
Should we hardcode to exclude the ground state?
| if seed is not None: | ||
| np.random.seed(seed) | ||
| valid = [] | ||
| temperature = 300 |
There was a problem hiding this comment.
The input parameter temp does not work, since temperature = 300 is hardcoded.
| pop.append(1 - np.sum(pop)) | ||
|
|
||
| if seed is not None: | ||
| np.random.seed(seed) |
There was a problem hiding this comment.
If the seed is set here, each wignerfunc call
q = np.random.uniform(0, 1) * 10.0 - 5.0
p = np.random.uniform(0, 1) * 10.0 - 5.0
wii give the same result.
3 participants
A follow-up to PR #493.