Serenity interface python code and examples

CHANGELOG.md

@@ -27,6 +27,7 @@ This changelog is effective from the 2025 releases.
 * Logging of job summaries to CSV logfile
 * Logging of AMS job error messages to stdout and logfile on job failure
 * Method `get_errormsg` enforced on the `Job` base class, with a default implementation
+* Added an interface to the Serenity program through methods such as `SerenityJob`, `SerenityResults` and `SerenitySettings`
 
 ### Changed
 * Functions for optional packages (e.g. RDKit, ASE) are available even when these packages are not installed, but will raise an `MissingOptionalPackageError` when called

doc/source/conf.py

@@ -312,6 +312,10 @@ def setup(app):
 .. |VASPJob| replace:: :class:`VASPJob<scm.plams.interfaces.thirdparty.vasp.VASPJob>`
 .. |VASPResults| replace:: :class:`VASPResults<scm.plams.interfaces.thirdparty.vasp.VASPResults>`
 
+.. |SerenityJob| replace:: :class:`~scm.plams.interfaces.serenity.SerenityJob`
+.. |SerenitySettings| replace:: :class:`~scm.plams.interfaces.serenity.SerenitySettings`
+.. |SerenityResults| replace:: :class:`~scm.plams.interfaces.serenity.SerenityResults`
+
 .. |CRSJob| replace:: :class:`CRSJob<scm.plams.interfaces.adfsuite.crs.CRSJob>`
 .. |CRSResults| replace:: :class:`CRSResults<scm.plams.interfaces.adfsuite.crs.CRSResults>`
 

doc/source/general.rst

@@ -113,6 +113,7 @@ Added
 * Support for AMS ``ChemicalSystem`` within |AMSJob| and |AMSResults|. |AMSJob| can accept a ``ChemicalSystem`` as an input system, and the methods :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_system`, :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_input_system` and :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_main_system` on |AMSResults| return a ``ChemicalSystem``. These provide the option to use a ``ChemicalSystem`` in place of a PLAMS ``Molecule``.
 * Support for work functions through :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_work_function_results` and :func:`~scm.plams.tools.plot.plot_work_function`
 * Added :meth:`~scm.plams.interfaces.molecule.packmol.packmol_around` method to pack around an existing molecule or pack molecules into a non-orthorhombic box.
+* Added an interface to the `Serenity program <https://github.com/qcserenity/serenity>`_ through classes such as |SerenityJob|, |SerenitySettings|, and |SerenityResults|. 
 
 Changed
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

doc/source/interfaces/serenity.rst

@@ -0,0 +1,89 @@
+Serenity
+-------------------------
+
+.. currentmodule:: scm.plams.interfaces.thirdparty.serenity
+
+Serenity is a quantum chemistry code originally developed in the group of Johannes Neugebauer at the University of Münster with a strong focus on quantum chemical subsystem/embedding methods. See the `serenity page <https://github.com/qcserenity/serenity>`_ for more details.
+
+PLAMS offers a simple Serenity interface and is capable of running Serenity calculations.
+The relevant classes are |SerenityJob|, |SerenitySettings|, and |SerenityResults|.
+
+It is also possible to run post-Hartree-Fock (HF) Serenity calculations starting from a ground-state Self-Consistent-Field calculation perfromed with ADF.
+
+
+Serenity calculation
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Preparing an instance of a |SerenityJob| follows general principles for |SingleJob|.
+Here is an example of a simple HF calculation::
+
+    from scm.plams.interfaces.thirdparty.serenity import SerenitySettings, SerenityJob
+
+    sersett = SerenitySettings()
+    # The path to the geometry file
+    sersett.input.system.gly_gly_gly.geometry = "mygeometry.xyz"
+    # Charge and spin
+    sersett.input.system.gly_gly_gly.charge = "0"
+    sersett.input.system.gly_gly_gly.spin = "0"
+    # The electronic structure method
+    sersett.input.system.gly_gly_gly.method = "HF"
+    sersett.input.system.gly_gly_gly.basis.label = "6-31G"
+    # Performs a SCF calculation 
+    sersett.input.task.SCF.act = "mymolecule"
+
+    serjob = SerenityJob(settings=sersett, name="Serenitycalc")
+    serjob.run()
+
+This input performs an HF/6-31G calculation on the chosen molecule.
+
+
+ADF interface
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Interfacing ADF with Serenity is done through a special file, saved by ADF, which contains information about the Hamiltonian matrix which can be used by Serenity for post-HF calculations.
+Here is a simple example::
+
+
+    from scm.plams import init, Molecule, Settings, AMSJob
+    from scm.plams.interfaces.thirdparty.serenity import SerenitySettings, SerenityJob
+
+    init(folder="PLAMS_CCSD_N2_DZ")
+
+    mol = Molecule("N2.xyz")
+    adfsett = Settings()
+    adfsett.input.ams.task = "SinglePoint"
+    adfsett.input.adf.IntegralsToFile = "SERENITY"
+    adfsett.input.adf.TotalEnergy = ""
+    adfsett.input.adf.basis.core = "None"
+    adfsett.input.adf.basis.type = "DZ"
+    adfsett.input.adf.basis.CreateOutput = "yes"
+    adfsett.input.adf.relativity = "Level=None"
+    adfsett.input.adf.symmetry = "NoSym"
+    adfsett.input.adf.xc.hartreefock = ""
+
+    adfjob = AMSJob(molecule=mol, settings=adfsett, name="ADF_N2_HF")
+    adfjob.run()
+    bonding_energy = adfjob.results.get_energy()
+    print(f"ADF bonding energy: {bonding_energy} hartree")
+
+
+    sersett = SerenitySettings()
+    sersett.input.system.N2.geometry = "N2.xyz"
+
+    sersett.input.task.CC.system = "N2"
+    sersett.input.task.CC.level = "CCSD"
+
+    serjob = SerenityJob(settings=sersett, name="Serenity_N2_CCSD")
+    serjob.run()
+    ccsd_correction = serjob.results.get_ccsd_energy_correction()
+    print(f"Serenity CCSD energy correction: {ccsd_correction} hartree")
+
+
+This script produces the CCSD/DZ energy of the nitrogen molecules using integrals computed by ADF with a Slater-type orbital basis.
+
+API
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. autoclass:: SerenityJob
+    :exclude-members: _result_type, _get_ready
+.. autoclass:: SerenityResults

doc/source/interfaces/thirdparty.rst

@@ -20,3 +20,4 @@ They offer a way of using PLAMS |Molecule| class with other libraries capable of
     orca
     raspa
     vasp
+    serenity

examples/SerenityAMS_DataTransfer_examples/C2H4.xyz

@@ -0,0 +1,8 @@
+6
+Ethene
+  H      1.1610      0.0661      1.0238
+  C      0.6579     -0.0045      0.0639
+  H      1.3352     -0.0830     -0.7815
+  C     -0.6579      0.0045     -0.0639
+  H     -1.3355      0.0830      0.7812
+  H     -1.1608     -0.0661     -1.0239

examples/SerenityAMS_DataTransfer_examples/C2H4_CCSD.py

@@ -0,0 +1,34 @@
+from scm.plams import init, Molecule, Settings, AMSJob
+from scm.plams.interfaces.thirdparty.serenity import SerenitySettings, SerenityJob
+
+init(folder="PLAMS_CCSD_C2H4_DZ")
+
+mol = Molecule("C2H4.xyz")
+adfsett = Settings()
+adfsett.input.ams.task = "SinglePoint"
+adfsett.input.adf.IntegralsToFile = "SERENITY"
+adfsett.input.adf.TotalEnergy = ""
+adfsett.input.adf.NumericalQuality = "VeryGood"
+adfsett.input.adf.basis.core = "None"
+adfsett.input.adf.basis.type = "DZ"
+adfsett.input.adf.basis.CreateOutput = "yes"
+adfsett.input.adf.relativity = "Level=None"
+adfsett.input.adf.symmetry = "NoSym"
+adfsett.input.adf.xc.hartreefock = ""
+
+adfjob = AMSJob(molecule=mol, settings=adfsett, name="ADF_C2H4_HF")
+adfjob.run()
+bonding_energy = adfjob.results.get_energy()
+print(f"ADF bonding energy: {bonding_energy} hartree")
+
+
+sersett = SerenitySettings()
+sersett.input.system.C2H4.geometry = "C2H4.xyz"
+
+sersett.input.task.CC.system = "C2H4"
+sersett.input.task.CC.level = "CCSD"
+
+serjob = SerenityJob(settings=sersett, name="Serenity_C2H4_CCSD")
+serjob.run()
+ccsd_correction = serjob.results.get_ccsd_energy_correction()
+print(f"Serenity CCSD energy correction: {ccsd_correction} hartree")

examples/SerenityAMS_DataTransfer_examples/HF.xyz

@@ -0,0 +1,4 @@
+2
+
+   F          0.91860        0.00000      500.00000
+   H          0.00000        0.00000      500.00000

examples/SerenityAMS_DataTransfer_examples/HF_CCSD_T.py

@@ -0,0 +1,34 @@
+from scm.plams import init, Molecule, Settings, AMSJob
+from scm.plams.interfaces.thirdparty.serenity import SerenitySettings, SerenityJob
+
+init(folder="PLAMS_CCSD_T_HF_TZP")
+
+mol = Molecule("HF.xyz")
+adfsett = Settings()
+adfsett.input.ams.task = "SinglePoint"
+adfsett.input.adf.IntegralsToFile = "SERENITY"
+adfsett.input.adf.TotalEnergy = ""
+adfsett.input.adf.NumericalQuality = "VeryGood"
+adfsett.input.adf.basis.core = "None"
+adfsett.input.adf.basis.type = "TZP"
+adfsett.input.adf.basis.CreateOutput = "yes"
+adfsett.input.adf.relativity = "Level=None"
+adfsett.input.adf.symmetry = "NoSym"
+adfsett.input.adf.xc.hartreefock = ""
+
+adfjob = AMSJob(molecule=mol, settings=adfsett, name="ADF_HF_HF")
+adfjob.run()
+bonding_energy = adfjob.results.get_energy()
+print(f"ADF bonding energy: {bonding_energy} hartree")
+
+
+sersett = SerenitySettings()
+sersett.input.system.HF.geometry = "HF.xyz"
+
+sersett.input.task.CC.system = "HF"
+sersett.input.task.CC.level = "CCSD(T)"
+
+serjob = SerenityJob(settings=sersett, name="Serenity_HF_CCSD_T")
+serjob.run()
+ccsd_correction = serjob.results.get_ccsd_energy_correction()
+print(f"Serenity CCSD energy correction: {ccsd_correction} hartree")

examples/SerenityAMS_DataTransfer_examples/N2.xyz

@@ -0,0 +1,4 @@
+2
+nitrogen
+  N         0.000000        0.000000        0.000000
+  N         1.097600        0.000000        0.000000

examples/SerenityAMS_DataTransfer_examples/N2_CCSD.py

@@ -0,0 +1,34 @@
+from scm.plams import init, Molecule, Settings, AMSJob
+from scm.plams.interfaces.thirdparty.serenity import SerenitySettings, SerenityJob
+
+init(folder="PLAMS_CCSD_N2_DZ")
+
+mol = Molecule("N2.xyz")
+adfsett = Settings()
+adfsett.input.ams.task = "SinglePoint"
+adfsett.input.adf.IntegralsToFile = "SERENITY"
+adfsett.input.adf.TotalEnergy = ""
+adfsett.input.adf.NumericalQuality = "VeryGood"
+adfsett.input.adf.basis.core = "None"
+adfsett.input.adf.basis.type = "DZ"
+adfsett.input.adf.basis.CreateOutput = "yes"
+adfsett.input.adf.relativity = "Level=None"
+adfsett.input.adf.symmetry = "NoSym"
+adfsett.input.adf.xc.hartreefock = ""
+
+adfjob = AMSJob(molecule=mol, settings=adfsett, name="ADF_N2_HF")
+adfjob.run()
+bonding_energy = adfjob.results.get_energy()
+print(f"ADF bonding energy: {bonding_energy} hartree")
+
+
+sersett = SerenitySettings()
+sersett.input.system.N2.geometry = "N2.xyz"
+
+sersett.input.task.CC.system = "N2"
+sersett.input.task.CC.level = "CCSD"
+
+serjob = SerenityJob(settings=sersett, name="Serenity_N2_CCSD")
+serjob.run()
+ccsd_correction = serjob.results.get_ccsd_energy_correction()
+print(f"Serenity CCSD energy correction: {ccsd_correction} hartree")

examples/SerenityCalculator/example.py

@@ -0,0 +1,30 @@
+import os
+from scm.plams import from_smiles, Settings, AMSJob, log
+
+molecule = from_smiles("[HH]")
+
+# first do a pesscan with pure AMS
+s = Settings()
+s.input.ams.Task = "PESScan"
+s.input.ams.PESScan.ScanCoordinate.nPoints = 5
+s.input.ams.PESScan.ScanCoordinate.Distance = "1 2 0.65 0.8"
+s.input.ForceField.Type = "UFF"
+psjob = AMSJob(settings=s, molecule=molecule, name="ams_pesscan")
+psjob.run()
+
+# then run a replay with the serenity calculator
+
+s = Settings()
+s.input.ams.Task = "Replay"
+s.input.ams.Replay.File = psjob.results.rkfpath()
+s.input.ASE.File = os.path.abspath("my_calculator.py")
+
+job = AMSJob(settings=s, molecule=molecule, name="serenity_replay")
+job.run()
+
+errormsg = job.get_errormsg()
+if errormsg:
+    log(errormsg)
+
+replayresults = job.results.get_pesscan_results()
+log(replayresults)

examples/SerenityCalculator/my_calculator.py

@@ -0,0 +1,14 @@
+from scm.plams.interfaces.thirdparty.serenity import SerenitySettings
+from scm.plams.interfaces.thirdparty.serenity_calculator import SerenityCalculator
+
+
+def get_calculator() -> SerenityCalculator:
+    s = SerenitySettings()
+    s.input.system.A.charge = "0"
+    s.input.system.A.spin = "0"
+    s.input.system.A.method = "HF"
+    s.input.system.A.basis.label = "6-31GS"
+    s.input.task.SCF.act = "A"
+
+    calc = SerenityCalculator(settings=s)
+    return calc

examples/SerenityInterface_examples/dft_res.py

@@ -0,0 +1,31 @@
+from scm.plams import init
+from scm.plams.interfaces.thirdparty.serenity import SerenitySettings, SerenityJob
+
+init(folder="test_riboflavin")
+
+# instead of giving the geometry key and the file as value you can also do:
+# Riboflavin = Molecule()
+# Riboflavin.add_atom(Atom(symbol='C', coords=(-2.335,1.498,5.390)))
+# Riboflavin.add_atom(Atom(symbol='C', coords=(-2.573,0.058,4.841)))
+# Riboflavin.add_atom(Atom(symbol='C', coords=(-3.154,0.011,3.393)))
+# ...
+# name of the molecule does not necessarily have to match the name of the system. You have to give the molecule to the job
+# serjob = SerenityJob(molecule=Riboflavin, settings=sersett, name='Serenity_Riboflavin')
+
+
+sersett = SerenitySettings()
+# The path to the geometry file
+sersett.input.system.Riboflavin.geometry = "riboflavin.xyz"
+# The electronic structure method
+sersett.input.system.Riboflavin.method = "DFT"
+sersett.input.system.Riboflavin.basis.label = "6-31GS"
+# Settings used for the DFT calculation.
+sersett.input.system.Riboflavin.dft.functional = "PBE"
+# Performs a SCF calculation for the given active (act) system
+sersett.input.task.SCF.act = "Riboflavin"
+# Every task in Serenity has a print-level which can be tuned to
+# adjust the amount of output generated.
+sersett.input.task.SCF.printLevel = "normal"
+
+serjob = SerenityJob(settings=sersett, name="Serenity_Riboflavin")
+serjob.run()

examples/SerenityInterface_examples/fde.py

@@ -0,0 +1,29 @@
+from scm.plams import init, Molecule, Atom
+from scm.plams.interfaces.thirdparty.serenity import SerenitySettings, SerenityJob
+
+init(folder="test")
+
+mol1 = Molecule()
+mol1.add_atom(Atom(symbol="O", coords=(0, 0, 0)))
+mol1.add_atom(Atom(symbol="H", coords=(0.758602, 0, 0.504284)))
+mol1.add_atom(Atom(symbol="H", coords=(0.260455, 0, -0.872893)))
+
+mol2 = Molecule()
+mol2.add_atom(Atom(symbol="O", coords=(3, 0.5, 0)))
+mol2.add_atom(Atom(symbol="H", coords=(3.758602, 0.5, 0.504284)))
+mol2.add_atom(Atom(symbol="H", coords=(3.260455, 0.5, -0.872893)))
+
+sersett = SerenitySettings()
+sersett.input.system.water1.method = "dft"
+sersett.input.system.water1.dft.functional = "pw91"
+
+sersett.input.system.water2.method = "dft"
+sersett.input.system.water2.dft.functional = "pw91"
+
+sersett.input.task.fde.act = "water1"
+sersett.input.task.fde.env = "water2"
+sersett.input.task.fde.emb.naddxcfunc = "pw91"
+sersett.input.task.fde.emb.naddkinfunc = "pw91k"
+# when using a dictionary and the molecule class to provide geometries you have no make sure that the given names match the systems that should use the respective geometry
+serjob = SerenityJob(molecule={"water1": mol1, "water2": mol2}, settings=sersett, name="water_dimer")
+serjob.run()

examples/SerenityInterface_examples/gly-gly-gly.xyz

@@ -0,0 +1,26 @@
+24
+gly-gly-gly.out	Energy:    -174.7655800
+N         -2.04230       -2.89390        2.81300
+H         -1.13560       -3.27750        2.63740
+C         -2.02210       -1.44800        2.54880
+H         -2.71180       -3.33290        2.21390
+C         -1.63580       -1.20010        1.07880
+H         -2.99220       -1.03750        2.73680
+H         -1.30570       -0.97820        3.18990
+O         -0.88210       -2.01130        0.48090
+N         -2.14750       -0.01590        0.37400
+H         -2.23950       -0.22170       -0.60020
+C         -1.21100        1.10310        0.55210
+C         -1.74710        2.34360       -0.18620
+H         -1.11250        1.32330        1.59460
+H         -0.25520        0.83540        0.15230
+O         -2.98730        2.50950       -0.31940
+N         -0.80990        3.33410       -0.73540
+H          0.01930        3.35400       -0.17680
+C         -1.43940        4.66250       -0.72550
+C         -0.45770        5.70010       -1.30080
+H         -2.32670        4.64120       -1.32320
+H         -1.69250        4.92970        0.27920
+O          0.82490        5.52860       -1.16310
+O         -0.91070        6.74840       -1.92470
+H         -0.29860        7.39520       -2.28340