SOC (1e and 2e cases)

README.md

@@ -8,6 +8,7 @@ Open Quantum Platform ([OpenQP](https://pubs.acs.org/doi/10.1021/acs.jctc.4c0111
 - **Flexible prototyping through a Python wrapper, PyOQP**
 - **Ground and Excited State Properties** by [MRSF-TDDFT](https://doi.org/10.1021/acs.jpclett.3c02296)
 - **Nonadiabatic Coupling** based on [TLF Technology](https://doi.org/10.1021/acs.jpclett.1c00932) using **MRSF-TDDFT**
+- **Spin-Orbit Coupling** by [**Relativistic** MRSF-TDDFT](https://doi.org/10.1021/acs.jctc.2c01036), including one-electron and mean-field two-electron SOC terms through `runtype=soc`
 - **New Exchange-Correlation Functionals** of [**DTCAM** series](https://doi.org/10.1021/acs.jctc.4c00640) for MRSF-TDDFT
 - **Ground State Properties** by HF and DFT theories
 - **Geometry Optimization, Transition State Search, and Conical Intersection Search** by SciPy and [DL-Find](https://github.com/digital-chemistry-laboratory/libdlfind)
@@ -23,7 +24,6 @@ Open Quantum Platform ([OpenQP](https://pubs.acs.org/doi/10.1021/acs.jctc.4c0111
 
 ### Upcoming Features
 - **Efficient electrostatic embedding QM/MM** by [ESPF QM/MM](https://doi.org/10.1063/5.0133646)
-- **Spin-Orbit Coupling** by [**Relativistic** MRSF-TDDFT](https://doi.org/10.1021/acs.jctc.2c01036)
 - **Ionization Potential/Electron Affinity** by [**EKT**-MRSF-TDDFT](https://doi.org/10.1021/acs.jpclett.1c02494)
 
 ### Quickstart

include/oqp.h

@@ -136,6 +136,8 @@ struct control_parameters {
     bool      basis_set_issue;
     double    conf_print_threshold;
     bool      rstctmo;
+    int64_t   scal_rel;
+    int64_t   soc_2e;  
     int64_t   converger_type;
     double    soscf_lvl_shift;
     int64_t   soscf_reset_mod;
@@ -225,4 +227,5 @@ void resp_charges(struct oqp_handle_t *inf);
 void mulliken(struct oqp_handle_t *inf);
 void mulliken_excited(struct oqp_handle_t *inf);
 void lowdin(struct oqp_handle_t *inf);
-
+void soc_mrsf(struct oqp_handle_t *inf);
+void dk_scalar(struct oqp_handle_t *inf);

pyoqp/README.md

@@ -74,6 +74,7 @@ Results, including log files and `test_report.txt`, will be stored in the curren
     runtype=energy
     system=''
     d4=False
+    soc_2e=1
 
 [guess]
     type=huckel
@@ -112,7 +113,6 @@ Results, including log files and `test_report.txt`, will be stored in the curren
     td_prop=False
     grad=0
     nac=''
-    soc=''
     export=False
 
 [optimize]
@@ -183,7 +183,7 @@ input section handle the basic information of molecular system
       grad       single-point energy and gradients
       hess       frequency calculation (numerical only)
       nac        non-adiabatic coupling (not available yet)
-      soc        spin-orbit coupling (not available yet)
+      soc        spin-orbit coupling for MRSF-TDDFT
       optimize   local minimum geometry optimization
       meci       minimum energy conical intersection optimization
       mep        minimum energy path calculation
@@ -207,6 +207,14 @@ input section handle the basic information of molecular system
       False      do not compute DFTD4 corrections according to functional (default)
       True       compute DFTD4 corrections for energy and gradients. some functional might not be supported
 
+- soc_2e // choose SOC terms for runtype=soc
+
+      0          one-electron SOC only
+      1          one-electron plus mean-field two-electron SOC (default)
+
+      runtype=soc requires method=tdhf, tdhf.type=mrsf, scf.type=rohf,
+      and scf.multiplicity=3.
+
 ### [guess]
 
 guess section handle the guess orbitals 

pyoqp/oqp/library/runfunc.py

@@ -87,7 +87,27 @@ def compute_nac(mol):
     NAC(mol).nac()
 
 def compute_soc(mol):
-    pass
+    sp = SinglePoint(mol)
+    ref_energy = sp.reference()          # SCF один раз
+
+    mol.data.set_tdhf_multiplicity(1)
+    mol.singlet_energies = sp.excitation(ref_energy)
+
+    mol.data['OQP::td_singlet_energies'] = mol.data['OQP::td_energies']
+    mol.data['OQP::td_bvec_mo_s'] = mol.data['OQP::td_bvec_mo'].copy()
+
+    mol.data.set_tdhf_multiplicity(3)
+    mol.triplet_energies = sp.excitation(ref_energy)
+
+    mol.data['OQP::td_triplet_energies'] = mol.data['OQP::td_energies']
+    mol.data['OQP::td_bvec_mo_t'] = mol.data['OQP::td_bvec_mo'].copy()
+
+    oqp.soc_mrsf(mol)
+
+    LastStep(mol).compute(mol) 
+
+#    mol.data['OQP::soc_']
+
 
 
 def compute_hess(mol):
@@ -143,11 +163,11 @@ def compute_properties(mol):
         pass
 
     # compute soc
-    soc_type = mol.config['properties']['soc']
-    if soc_type:
-        pass
-    else:
-        pass
+#    soc_type = mol.config['properties']['soc']
+#    if soc_type:
+#        pass
+#    else:
+#        pass
 
 def compute_data(mol):
     # compute reference energy
@@ -167,11 +187,11 @@ def compute_data(mol):
         pass
 
     # compute soc
-    soc_type = mol.config['properties']['soc']
-    if soc_type:
-        pass
-    else:
-        pass
+#    soc_type = mol.config['properties']['soc']
+#    if soc_type:
+#        pass
+#    else:
+#        pass
 
 
 def get_optimizer(mol):

pyoqp/oqp/molecule/molecule.py

@@ -59,6 +59,11 @@ def __init__(self, project_name, input_file, log,
             'OQP::td_abxc', 'OQP::td_bvec_mo', 'OQP::td_mrsf_density', 'OQP::td_energies',
             'OQP::td_states_overlap',
             'OQP::dc_matrix', 'OQP::nac_matrix',
+            'OQP::td_singlet_energies', 'OQP::td_triplet_energies',
+            'OQP::td_bvec_mo_s', 'OQP::td_bvec_mo_t',
+            'OQP::soc_eval',
+            'OQP::soc_evec_re', 'OQP::soc_evec_im',
+            'OQP::soc_hsoc_re', 'OQP::soc_hsoc_im',
         ]
         self.skip_tag = {"rhf": ['OQP::DM_B', 'OQP::FOCK_B', 'OQP::E_MO_B', 'OQP::VEC_MO_B'],
                          "rohf": [],

pyoqp/oqp/molecule/oqpdata.py

@@ -55,6 +55,7 @@ def path(strng):
         'system': {'type': str, 'default': ''},
         'system2': {'type': str, 'default': ''},
         'd4': {'type': bool, 'default': 'False'},
+        'soc_2e': {'type': int,  'default': '1'},    
     },
     'guess': {
         'type': {'type': string, 'default': 'huckel'},
@@ -95,6 +96,7 @@ def path(strng):
         'save_molden': {'type': bool, 'default': 'True'},
         'rstctmo': {'type': bool, 'default': 'False'},
         'converger_type': {'type': string, 'default': 'diis'},
+        'scal_rel': {'type': int, 'default': '0'},
         'soscf_reset_mod': {'type': int, 'default': '0'},
         'soscf_mode': {'type': int, 'default': '0'},
         'soscf_lvl_shift': {'type': float, 'default': '0'},
@@ -153,7 +155,6 @@ def path(strng):
         'td_prop': {'type': bool, 'default': 'False'},
         'grad': {'type': iarray, 'default': '0'},
         'nac': {'type': str, 'default': ''},
-        'soc': {'type': str, 'default': ''},
         'export': {'type': bool, 'default': 'False'},
         'title': {'type': str, 'default': ''},
         'back_door': {'type': bool, 'default': False}
@@ -255,6 +256,7 @@ class OQPData:
             "functional": "set_dft_functional",
             "system": "set_system",
             "system2": "set_system2",
+            "soc_2e":     "set_soc_2e",
         },
         "guess": {
         },
@@ -282,6 +284,7 @@ class OQPData:
             "incremental": "set_scf_incremental",
             "active_basis": "set_scf_active_basis",
             "rstctmo": "set_scf_rstctmo",
+            "scal_rel": "set_scf_scal_rel",
             "converger_type": "set_scf_converger_type",
             "soscf_reset_mod": "set_scf_soscf_reset_mod",
             "soscf_mode": "set_scf_soscf_mode",
@@ -563,6 +566,10 @@ def set_scf_conv(self, conv):
         """Set SCF convergence threshold"""
         self._data.control.conv = conv
 
+    def set_scf_scal_rel(self, scal_rel):
+        """Set SCF convergence threshold"""
+        self._data.control.scal_rel = scal_rel
+
     def set_scf_incremental(self, flag):
         """Set incremental Fock matrix build"""
         self._data.control.scf_incremental = 1 if flag else 0
@@ -855,6 +862,9 @@ def set_system2(self, system):
             num_atoms, x, y, z, q, mass = read_system(system)
             self.mol2 = np.array(x + y + z).reshape((3, num_atoms)).T.reshape(-1)
 
+    def set_soc_2e(self, soc_2e):
+        self._data.control.soc_2e = soc_2e
+
     def parse_section(self, config, section):
         cfg_input = config[section]
         for key in cfg_input.keys():

pyoqp/oqp/utils/input_checker.py

@@ -12,9 +12,9 @@
 
 SUPPORTED_RUNTYPES = {
     "energy", "grad", "hess", "nac", "nacme", "bp", "optimize",
-    "meci", "mecp", "mep", "ts", "irc", "prop", "data",
+    "meci", "mecp", "mep", "ts", "irc", "prop", "data", "soc",
 }
-NOT_AVAILABLE_RUNTYPES = {"soc", "neb"}
+NOT_AVAILABLE_RUNTYPES = {"neb"}
 ALL_RUNTYPES = SUPPORTED_RUNTYPES | NOT_AVAILABLE_RUNTYPES
 
 METHODS = {"hf", "tdhf"}
@@ -36,7 +36,7 @@
 INIT_SCF_TYPES = {"no", "rhf", "uhf", "rohf", "rks", "uks", "roks"}
 
 WIKI_HELP = {
-    "input.runtype": "Use energy, grad, hess, nac, nacme, optimize, meci, mecp, mep, ts, prop, or data. soc and neb are recognized but not implemented yet.",
+    "input.runtype": "Use energy, grad, hess, nac, nacme, soc, optimize, meci, mecp, mep, ts, prop, or data. neb is recognized but not implemented yet.",
     "input.method": "Use method=hf for HF/DFT and method=tdhf for TDHF/TDDFT/SF/MRSF runs.",
     "input.system": "Set system to an XYZ file path or inline coordinates with one atom per indented line.",
     "input.basis": "Set basis to a basis name, a comma-separated per-atom list, or library with tagged atoms and [input] library mappings.",
@@ -738,6 +738,9 @@ def _check_runtype(config: dict[str, Any], report: CheckReport) -> None:
     if runtype == "nacme":
         _check_nacme(config, report)
 
+    if runtype == "soc":       
+        _check_soc(config, report)
+
     if runtype == "hess":
         _check_hess(config, report)
 
@@ -962,6 +965,51 @@ def _check_dlfind(config: dict[str, Any], report: CheckReport) -> None:
                 action="Use icoord 0-4 for TS optimization.",
             )
 
+def _check_soc(config: dict[str, Any], report: CheckReport) -> None:
+    method = _as_lower(_get(config, "input", "method", "hf"))
+    td_type = _as_lower(_get(config, "tdhf", "type", "rpa"))
+    scf_type = _as_lower(_get(config, "scf", "type", "rhf"))
+    scf_mult = _get(config, "scf", "multiplicity", 1)
+
+    if method != "tdhf":
+        report.add(
+            "ERROR",
+            "input.method",
+            "SOC requires method=tdhf.",
+            value=method,
+            expected="tdhf",
+            action="Set [input] method=tdhf.",
+        )
+
+    if td_type != "mrsf":
+        report.add(
+            "ERROR",
+            "tdhf.type",
+            "SOC requires tdhf.type=mrsf.",
+            value=td_type,
+            expected="mrsf",
+            action="Set [tdhf] type=mrsf.",
+        )
+
+    if scf_type != "rohf":
+        report.add(
+            "ERROR",
+            "scf.type",
+            "SOC requires an ROHF reference.",
+            value=scf_type,
+            expected="rohf",
+            action="Set [scf] type=rohf.",
+        )
+
+    if scf_mult != 3:
+        report.add(
+            "ERROR",
+            "scf.multiplicity",
+            "SOC requires a triplet ROHF reference (multiplicity=3).",
+            value=scf_mult,
+            expected="3",
+            action="Set [scf] multiplicity=3.",
+        )
 
 def _check_hess(config: dict[str, Any], report: CheckReport) -> None:
     method = _as_lower(_get(config, "input", "method", "hf"))