Operators

README.md

@@ -3,15 +3,30 @@ ResponseFun
 
 Fun with Response Functions in the Algebraic Diagrammatic Construction Framework.
 
+
+### Installation
+
+You can install `responsefun` with Conda using the conda-forge channel:
+
+```bash
+conda install responsefun -c conda-forge
+```
+
+
 ### Development
 
 A suitable `conda` environment for development can be created from the `ci_env.yml` file.
+After cloning the repository and navigating to the responsefun directory, the `responsefun` python package can be installed with the following command:
+`pip install -e .`
 
 Tests can be run with `pytest responsefun`, or `pytest --pyargs responsefun` if the package is installed.
-To exclude the slowest test, run `pytest -k "not slow" responsefun`.
+To exclude the slowest test, run `pytest -m "not slow" responsefun`.
 
 Code style is enforced through `black` (formatting), `isort` (sorting import statements), and `ruff` (linting).
 
+### Citation
+
+If you use `responsefun`, please cite [our article in JCTC.](https://doi.org/10.1021/acs.jctc.3c00456)
 ### Copyright
 
 Copyright (c) 2023, The `responsefun` Developers

ci_env.yml

@@ -3,10 +3,10 @@ channels:
   - conda-forge
 dependencies:
   - libblas =*=*mkl
-  - adcc >=0.15.16,<0.16.0
+  - adcc >=0.16.0
   - respondo >=0.0.5
   - numpy >=1.14
   - isort
   - ruff
   - pip:
-    - pyscf
+    - pyscf

examples/alpha.py

@@ -6,7 +6,7 @@
 from pyscf import gto, scf
 
 from responsefun import evaluate_property_isr, TransitionMoment
-from responsefun.symbols_and_labels import O, gamma, n, op_a, op_b, w, w_n
+from responsefun.symbols_and_labels import O, gamma, n, mu_a, mu_b, w, w_n
 
 # run SCF in PySCF
 mol = gto.M(
@@ -26,8 +26,8 @@
 
 # define symbolic SOS expression
 alpha_sos_expr = (
-        TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, O) / (w_n - w - 1j*gamma)
-        + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_a, O) / (w_n + w + 1j*gamma)
+        TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w - 1j*gamma)
+        + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, O) / (w_n + w + 1j*gamma)
 )
 # compute polarizability
 alpha_tens = evaluate_property_isr(

examples/beta_complex.py

@@ -9,9 +9,9 @@
     O,
     gamma,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     p,
     w_1,
     w_2,
@@ -38,19 +38,16 @@
 
 # compute the complex beta tensor
 beta_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n - w_o - 1j*gamma) * (w_p - w_2 - 1j*gamma))
-    + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_a, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n + w_1 + 1j*gamma) * (w_p - w_2 - 1j*gamma))
-    + TransitionMoment(O, op_b, n) * TransitionMoment(n, op_c, p, shifted=True)
-    * TransitionMoment(p, op_a, O) / ((w_n + w_1 + 1j*gamma) * (w_p + w_o + 1j*gamma))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n - w_o - 1j*gamma) * (w_p - w_2 - 1j*gamma))
+    + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_a, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n + w_1 + 1j*gamma) * (w_p - w_2 - 1j*gamma))
+    + TransitionMoment(O, mu_b, n) * TransitionMoment(n, mu_c, p, shifted=True)
+    * TransitionMoment(p, mu_a, O) / ((w_n + w_1 + 1j*gamma) * (w_p + w_o + 1j*gamma))
 )
-# the minus sign is needed, because the negative charge is not yet included
-# in the operator definitions
-# TODO: remove minus after adc-connect/adcc#190 is merged
-beta_tens = -1.0 * evaluate_property_isr(
+beta_tens = evaluate_property_isr(
     state, beta_term, [n, p],
-    perm_pairs=[(op_b, w_1), (op_c, w_2)], excluded_states=O,
+    perm_pairs=[(mu_b, w_1), (mu_c, w_2)], excluded_states=O,
     freqs_in=[(w_1, 0.5), (w_2, 0.5)], freqs_out=(w_o, w_1+w_2),
     damping=0.01, conv_tol=1e-4,
 )

examples/beta_sos.py

@@ -6,9 +6,9 @@
 from responsefun.symbols_and_labels import (
     O,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     p,
     w_1,
     w_2,
@@ -18,16 +18,16 @@
 )
 
 beta_sos_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n - w_o) * (w_p - w_2))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n - w_o) * (w_p - w_2))
 )
 
 beta_sos = SumOverStates(
     beta_sos_term,  # first SOS term
     [n, p],  # indices of summation
     freqs_in=[w_1, w_2],  # frequencies of incident photons
     freqs_out=w_o,  # frequency of resulting photon
-    perm_pairs=[(op_a, -w_o), (op_b, w_1), (op_c, w_2)],  # tuples to be permuted
+    perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)],  # tuples to be permuted
     excluded_states=O  # states excluded from the summations
 )
 

examples/esp.py

@@ -6,7 +6,7 @@
 from pyscf import gto, scf
 
 from responsefun import evaluate_property_isr, TransitionMoment
-from responsefun.symbols_and_labels import f, gamma, n, op_a, op_b, w, w_f, w_n
+from responsefun.symbols_and_labels import f, gamma, n, mu_a, mu_b, w, w_f, w_n
 
 # run SCF in PySCF
 mol = gto.M(
@@ -27,8 +27,8 @@
 
 # compute esp tensor
 sos_expr = (
-    TransitionMoment(f, op_a, n) * TransitionMoment(n, op_b, f) / (w_n - w_f - w - 1j*gamma)
-    + TransitionMoment(f, op_b, n) * TransitionMoment(n, op_a, f) / (w_n - w_f + w + 1j*gamma)
+    TransitionMoment(f, mu_a, n) * TransitionMoment(n, mu_b, f) / (w_n - w_f - w - 1j*gamma)
+    + TransitionMoment(f, mu_b, n) * TransitionMoment(n, mu_a, f) / (w_n - w_f + w + 1j*gamma)
 )
 tens = evaluate_property_isr(
     state,  # ExcitedStates object returned by the adcc calculation

examples/mcd.py

@@ -8,7 +8,7 @@
 
 from responsefun import evaluate_property_isr, TransitionMoment
 from responsefun.misc import epsilon
-from responsefun.symbols_and_labels import O, j, k, op_a, op_b, opm_c, w_j, w_k
+from responsefun.symbols_and_labels import O, j, k, mu_a, mu_b, m_c, w_j, w_k
 
 # run SCF in PySCF
 mol = gto.M(
@@ -28,11 +28,11 @@
 
 # define symbolic SOS expressions
 mcd_sos_expr1 = (
-        TransitionMoment(O, opm_c, k) * TransitionMoment(k, op_b, j, shifted=True)
-        * TransitionMoment(j, op_a, O) / w_k
+        TransitionMoment(O, m_c, k) * TransitionMoment(k, mu_b, j, shifted=True)
+        * TransitionMoment(j, mu_a, O) / w_k
 )
 mcd_sos_expr2 = (
-        TransitionMoment(O, op_b, k) * TransitionMoment(k, opm_c, j) * TransitionMoment(j, op_a, O)
+        TransitionMoment(O, mu_b, k) * TransitionMoment(k, m_c, j) * TransitionMoment(j, mu_a, O)
         / (w_k - w_j)
 )
 # compute MCD B term for the first excited state
@@ -48,8 +48,5 @@
     conv_tol=1e-4,
 )
 
-# the minus sign is needed, because the negative charge is not yet included
-# in the operator definitions
-# TODO: remove minus after adc-connect/adcc#190 is merged
-bterm = -1.0 * np.einsum("abc,abc->", epsilon, mcd_tens1 + mcd_tens2)
+bterm = np.einsum("abc,abc->", epsilon, mcd_tens1 + mcd_tens2)
 print(f"The MCD Faraday B term for excited state {final_state} is {bterm:.2f} (a.u.).")

examples/rixs.py

@@ -12,8 +12,8 @@
     f,
     gamma,
     n,
-    op_a,
-    op_b,
+    mu_a,
+    mu_b,
     w,
     w_f,
     w_n,
@@ -38,7 +38,7 @@
 state = adcc.adc2(scfres, n_singlets=5)
 
 # compute RIXS tensor within the rotating-wave approximation
-rixs_sos_rwa = TransitionMoment(f, op_a, n) * TransitionMoment(n, op_b, O) / (w_n - w - 1j * gamma)
+rixs_sos_rwa = TransitionMoment(f, mu_a, n) * TransitionMoment(n, mu_b, O) / (w_n - w - 1j * gamma)
 rixs_rwa = evaluate_property_isr(
     state,
     rixs_sos_rwa,
@@ -55,7 +55,7 @@
     state,
     rixs_sos_rwa,
     [n],
-    perm_pairs=[(op_a, w + 1j * gamma), (op_b, -w_prime - 1j * gamma)],
+    perm_pairs=[(mu_a, w + 1j * gamma), (mu_b, -w_prime - 1j * gamma)],
     freqs_in=(w, ev2au(534.74)),
     freqs_out=(w_prime, w-w_f),
     damping=ev2au(0.124),

examples/second_order_nlo.py

@@ -11,9 +11,9 @@
 from responsefun.symbols_and_labels import (
     O,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     p,
     w_1,
     w_2,
@@ -56,20 +56,17 @@ def compute_beta_parallel(beta_tens, dip_mom):
 
 # compute the first hyperpolarizability tensor
 beta_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, p, shifted=True)
-    * TransitionMoment(p, op_c, O) / ((w_n - w_o) * (w_p - w_2))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, p, shifted=True)
+    * TransitionMoment(p, mu_c, O) / ((w_n - w_o) * (w_p - w_2))
 )
 processes = {
     "static": (0.0, 0.0), "OR": (w_ruby, -w_ruby),
     "EOPE": (w_ruby, 0.0), "SHG": (w_ruby, w_ruby)
 }
 for process, freqs in processes.items():
-    # the minus sign is needed, because the negative charge is not yet included
-    # in the operator definitions
-    # TODO: remove minus after adc-connect/adcc#190 is merged
-    beta_tens = -1.0 * evaluate_property_isr(
+    beta_tens = evaluate_property_isr(
         state, beta_term, [n, p],
-        perm_pairs=[(op_a, -w_o), (op_b, w_1), (op_c, w_2)],
+        perm_pairs=[(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2)],
         freqs_in=[(w_1, freqs[0]), (w_2, freqs[1])],
         freqs_out=(w_o, w_1+w_2),
         excluded_states=O,

examples/third_order_nlo.py

@@ -12,10 +12,10 @@
     O,
     m,
     n,
-    op_a,
-    op_b,
-    op_c,
-    op_d,
+    mu_a,
+    mu_b,
+    mu_c,
+    mu_d,
     p,
     w_1,
     w_2,
@@ -55,16 +55,16 @@ def compute_gamma_average(gamma_tens):
 state = adcc.adc2(scfres, n_singlets=5)
 # compute the second hyperpolarizability tensor
 gamma_term_I = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, m, shifted=True)
-    * TransitionMoment(m, op_c, p, shifted=True) * TransitionMoment(p, op_d, O)
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, m, shifted=True)
+    * TransitionMoment(m, mu_c, p, shifted=True) * TransitionMoment(p, mu_d, O)
     / ((w_n - w_o) * (w_m - w_2 - w_3) * (w_p - w_3))
 )
 gamma_term_II = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, O)
-    * TransitionMoment(O, op_c, m) * TransitionMoment(m, op_d, O)
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, O)
+    * TransitionMoment(O, mu_c, m) * TransitionMoment(m, mu_d, O)
     / ((w_n - w_o) * (w_m - w_3) * (w_m + w_2))
 )
-perm_pairs = [(op_a, -w_o), (op_b, w_1), (op_c, w_2), (op_d, w_3)]
+perm_pairs = [(mu_a, -w_o), (mu_b, w_1), (mu_c, w_2), (mu_d, w_3)]
 processes = {
     "static": (0.0, 0.0, 0.0), "dcOR": (w_ruby, -w_ruby, 0.0),
     "EOKE": (w_ruby, 0.0, 0.0), "IDRI": (w_ruby, -w_ruby, w_ruby),

examples/threepa.py

@@ -11,9 +11,9 @@
     f,
     m,
     n,
-    op_a,
-    op_b,
-    op_c,
+    mu_a,
+    mu_b,
+    mu_c,
     w_1,
     w_2,
     w_3,
@@ -45,18 +45,15 @@ def threepa_average(tens):
 print(state.describe())
 
 threepa_term = (
-    TransitionMoment(O, op_a, n) * TransitionMoment(n, op_b, m)
-    * TransitionMoment(m, op_c, f) / ((w_n - w_1) * (w_m - w_1 - w_2))
+    TransitionMoment(O, mu_a, n) * TransitionMoment(n, mu_b, m)
+    * TransitionMoment(m, mu_c, f) / ((w_n - w_1) * (w_m - w_1 - w_2))
 )
 
 for es in range(5):
     print(f"===== State {es} ===== ")
-    # the minus sign is needed, because the negative charge is not yet included
-    # in the operator definitions
-    # TODO: remove minus after adc-connect/adcc#190 is merged
-    threepa_tens = -1.0 * evaluate_property_isr(
+    threepa_tens = evaluate_property_isr(
         state, threepa_term, [n, m],
-        perm_pairs=[(op_a, w_1), (op_b, w_2), (op_c, w_3)],
+        perm_pairs=[(mu_a, w_1), (mu_b, w_2), (mu_c, w_3)],
         freqs_in=[(w_1, w_f/3), (w_2, w_f/3), (w_3, w_f/3)],
         excited_state=es, conv_tol=1e-5,
     )

examples/tpa.py

@@ -10,8 +10,8 @@
     O,
     f,
     n,
-    op_a,
-    op_b,
+    mu_a,
+    mu_b,
     w_1,
     w_2,
     w_n,
@@ -42,14 +42,14 @@ def tpa_average(tens):
 print(state.describe())
 
 tpa_term = (
-    TransitionMoment(f, op_b, n) * TransitionMoment(n, op_a, O) / (w_n - w_1)
+    TransitionMoment(f, mu_b, n) * TransitionMoment(n, mu_a, O) / (w_n - w_1)
 )
 
 for es in range(1):
     print(f"===== State {es} ===== ")
     tpa_tens = evaluate_property_isr(
         state, tpa_term, [n],
-        perm_pairs=[(op_a, w_1), (op_b, w_2)],
+        perm_pairs=[(mu_a, w_1), (mu_b, w_2)],
         freqs_in=[(w_1, w_f/2), (w_2, w_f/2)],
         excited_state=es, conv_tol=1e-4,
     )

pyproject.toml

@@ -21,10 +21,10 @@ classifiers = [
 requires-python = ">=3.7"
 # Declare any run-time dependencies that should be installed with the package.
 dependencies = [
-    "adcc>=0.15.16,<0.16.0",
+    "adcc>=0.16.0",
     "respondo>=0.0.5",
     "numpy>=1.14",
-    "sympy",
+    "sympy<=1.13",
     "tqdm",
     "anytree",
 ]