Write and Read only nonzero Hamiltonian elements

include/hamiltonianfilereader.hpp

@@ -4,6 +4,7 @@
 #include <vector>
 #include <assert.h>
 #include <iostream> 
+#include <fstream>
 #include <controlbasis.hpp>
 
 #pragma once
@@ -20,8 +21,10 @@ class HamiltonianFileReader{
 
     LindbladType lindbladtype; ///< Type of solver (Lindblad vs Schroedinger)
     int dim_rho; ///< Dimension of the Hilbert space (N)
-    std::string hamiltonian_file; ///< Filename for Hamiltonian data ('none' if not used)
+    std::string hamiltonian_file_Hsys; ///< Filename for system Hamiltonian data ('none' if not used)
+    std::string hamiltonian_file_Hc; ///< Filename for control Hamiltonian data ('none' if not used)
     int mpirank_world; ///< Rank of global MPI communicator
+    int mpisize_world; ///< Rank of global MPI communicator
     bool quietmode; ///< Flag for quiet mode operation
 
 	public:
@@ -30,32 +33,32 @@ class HamiltonianFileReader{
     /**
      * @brief Constructor with Hamiltonian file specification.
      *
-     * @param hamiltonian_file_ Path to file containing Hamiltonian data
+     * @param hamiltonian_file_Hsys Path to file containing system Hamiltonian data
+     * @param hamiltonian_file_Hc Path to file containing control Hamiltonian data
      * @param lindbladtype_ Type of solver (Lindblad or Schroedinger)
      * @param dim_rho_ Dimension of the Hilbert space
      * @param quietmode_ Flag for quiet operation
      */
-    HamiltonianFileReader(std::string hamiltonian_file_, LindbladType lindbladtype_, int dim_rho_, bool quietmode_);
+    HamiltonianFileReader(std::string hamiltonian_file_Hsys, std::string hamiltonian_file_Hc, LindbladType lindbladtype_, int dim_rho_, bool quietmode_);
 
     ~HamiltonianFileReader();
 
   /**
    * @brief Reads the constant system Hamiltonian from file.
    *
-   * @param[out] Bd Reference to matrix that stores the imaginary part of the system matrix for (-i*Hsys). Must be allocated. 
    * @param[out] Ad Reference to matrix that stores the real part of the system matrix for (-i*Hsys). Must be allocated. 
+   * @param[out] Bd Reference to matrix that stores the imaginary part of the system matrix for (-i*Hsys). Must be allocated. 
    */
-  void receiveHsys(Mat& Bd, Mat& Ad);
+  void receiveHsys(Mat& Ad, Mat& Bd);
 
   /**
    * @brief Receives real and imaginary control operators from file.
    *
    * Reads control Hamiltonian matrices for each oscillator from the
    * Hamiltonian file and stores them in the provided matrix vectors.
    *
-   * @param noscillators Number of oscillators in the system
    * @param Ac_vec Reference to vector of matrices storing real parts of control system matrices. One per oscillator.
    * @param Bc_vec Reference to vector of matrices storing imaginary parts of control matrices. One per oscillator.
    */
-  void receiveHc(int noscillators, std::vector<Mat>& Ac_vec, std::vector<Mat>& Bc_vec);
+  void receiveHc(std::vector<Mat>& Ac_vec, std::vector<Mat>& Bc_vec);
 };

include/mastereq.hpp

@@ -115,7 +115,8 @@ class MasterEq{
     IS isu, isv; ///< Vector strides for accessing real and imaginary parts u=Re(x), v=Im(x)
     Vec aux; ///< Auxiliary vector for computations
     bool quietmode; ///< Flag for quiet mode operation
-    std::string hamiltonian_file; ///< Filename if a custom Hamiltonian is read from file ('none' if standard Hamiltonian is used)
+    std::string hamiltonian_file_Hsys; ///< Filename if a custom system Hamiltonian is read from file ('none' if standard Hamiltonian is used)
+    std::string hamiltonian_file_Hc; ///< Filename if a custom control Hamiltonians are read from file ('none' if standard Hamiltonian is used)
 
   public:
     std::vector<int> nlevels; ///< Number of levels per oscillator
@@ -137,10 +138,11 @@ class MasterEq{
      * @param eta_ Frequency differences for rotating frame
      * @param lindbladtype_ Type of Lindblad operators to include
      * @param usematfree_ Flag to use matrix-free solver
-     * @param hamiltonian_file Filename for Hamiltonian data
+     * @param hamiltonian_file_Hsys Filename for system Hamiltonian data
+     * @param hamiltonian_file_Hc Filename for control Hamiltonian data
      * @param quietmode Flag for quiet operation (default: false)
      */
-    MasterEq(const std::vector<int>& nlevels, const std::vector<int>& nessential, Oscillator** oscil_vec_, const std::vector<double>& crosskerr_, const std::vector<double>& Jkl_, const std::vector<double>& eta_, LindbladType lindbladtype_, bool usematfree_, const std::string& hamiltonian_file, bool quietmode=false);
+    MasterEq(const std::vector<int>& nlevels, const std::vector<int>& nessential, Oscillator** oscil_vec_, const std::vector<double>& crosskerr_, const std::vector<double>& Jkl_, const std::vector<double>& eta_, LindbladType lindbladtype_, bool usematfree_, const std::string& hamiltonian_file_Hsys, const std::string& hamiltonian_file_Hc, bool quietmode=false);
 
     ~MasterEq();
 

quandary.py

@@ -86,7 +86,8 @@ class Quandary:
     -------------------
     _ninit                : int  # number of initial conditions that are propagated
     _lindblad_solver      : bool # Flag to determine whether lindblad solver vs schroedinger solver 
-    _hamiltonian_filename : str 
+    _hamiltonian_filename_Hsys : str 
+    _hamiltonian_filename_Hc : str 
     _gatefilename         : str 
     _initstatefilename    : str 
     _initialstate         : List[complex] = field(default_factory=list)
@@ -159,7 +160,8 @@ class Quandary:
     # Internal configuration. Should not be changed by user.
     _ninit                : int         = -1
     _lindblad_solver      : bool        = False
-    _hamiltonian_filename : str         = ""
+    _hamiltonian_filename_Hsys : str         = ""
+    _hamiltonian_filename_Hc : str         = ""
     _gatefilename         : str         = ""
     _initstatefilename    : str         = ""
     _initialstate         : List[complex] = field(default_factory=list)
@@ -205,6 +207,7 @@ def __post_init__(self):
             self.initctrl_MHz = [10.0 for _ in range(len(self.Ne))]
         if len(self.Hsys) > 0 and not self.standardmodel: # User-provided Hamiltonian operators 
             self.standardmodel=False   
+            self.usematfree=False
         else: # Using standard Hamiltonian model. Set it up only if needed for computing dT or the carrier wave frequencies later
             self.standardmodel=True
         if len(self.targetstate) > 0:
@@ -557,35 +560,28 @@ def __dump(self, *, pcof0=[], runtype="simulation", datadir="./run_dir"):
         # If not standard Hamiltonian model, write provided Hamiltonians to a file
         if not self.standardmodel:
             # Write non-standard Hamiltonians to file  
-            self._hamiltonian_filename= "hamiltonian.dat"
-            with open(os.path.join(datadir, self._hamiltonian_filename), "w", newline='\n') as f:
-                f.write("# Hsys_real \n")
-                Hsyslist = list(np.array(self.Hsys.real).flatten(order='F'))
-                for value in Hsyslist:
-                    f.write("{:20.13e}\n".format(value))
-                f.write("# Hsys_imag\n")
-                Hsyslist = list(np.array(self.Hsys.imag).flatten(order='F'))
-                for value in Hsyslist:
-                    f.write("{:20.13e}\n".format(value))
-
-            # Write control Hamiltonians to file, if given (append to file)
-            for iosc in range(len(self.Ne)):
-                # Real part, if given
-                if len(self.Hc_re)>iosc and len(self.Hc_re[iosc])>0:
-                    with open(os.path.join(datadir, self._hamiltonian_filename), "a", newline='\n') as f:
-                        Hcrelist = list(np.array(self.Hc_re[iosc]).flatten(order='F'))
-                        f.write("# Oscillator {:d} Hc_real \n".format(iosc))
-                        for value in Hcrelist:
-                            f.write("{:20.13e}\n".format(value))
-                # Imaginary part, if given
-                if len(self.Hc_im)>iosc and len(self.Hc_im[iosc])>0:
-                    with open(os.path.join(datadir, self._hamiltonian_filename), "a", newline='\n') as f:
-                        Hcimlist = list(np.array(self.Hc_im[iosc]).flatten(order='F'))
-                        f.write("# Oscillator {:d} Hc_imag \n".format(iosc))
-                        for value in Hcimlist:
-                            f.write("{:20.13e}\n".format(value))
+            # Write system Hamiltonian (complex)
+            self._hamiltonian_filename_Hsys= "hamiltonian_Hsys.dat"
+            H = self.Hsys
+            with open(os.path.join(datadir, self._hamiltonian_filename_Hsys), "w", newline='\n') as f:
+                f.write("# row col Hsys_real Hsys_imag \n")
+                nz = np.nonzero(H)
+                for i, j in zip(*nz):
+                    v = H[i, j]
+                    f.write(f"{i} {j} {v.real:.13e} {v.imag:.13e}\n")
+            # Write control Hamiltonians, if given
+            if len(self.Hc_re)>0 or len(self.Hc_im)>0:
+                self._hamiltonian_filename_Hc = "hamiltonian_Hc.dat"
+                with open(os.path.join(datadir, self._hamiltonian_filename_Hc), "w", newline='\n') as f:
+                    for iosc, (Hc_re, Hc_im) in enumerate(zip(self.Hc_re, self.Hc_im)):
+                        Hc = np.array(Hc_re) + 1j * np.array(Hc_im)
+                        nz = np.nonzero(Hc)
+                        f.write(f"# oscillator row col Hc_real Hc_imag \n")
+                        for i, j in zip(*nz):
+                            v = Hc[i, j]
+                            f.write(f"{iosc} {i} {j} {v.real:.13e} {v.imag:.13e}\n")
             if self.verbose:
-                print("Hamiltonian operators written to ", os.path.join(datadir, self._hamiltonian_filename))
+                print("Hamiltonian operators written to ", os.path.join(datadir, self._hamiltonian_filename_Hsys), os.path.join(datadir, self._hamiltonian_filename_Hc))
 
         # Initializing the control parameter vector 'pcof0'
         # 1. If the initial parameter vector (list) is given with the 'pcof0' argument, the list will be dumped to a file with name self.pcof0_filename := "pcof0.dat". 
@@ -712,7 +708,10 @@ def __dump(self, *, pcof0=[], runtype="simulation", datadir="./run_dir"):
         mystring += "linearsolver_type = gmres\n"
         mystring += "linearsolver_maxiter = 20\n"
         if not self.standardmodel:
-            mystring += "hamiltonian_file= "+str(self._hamiltonian_filename)+"\n"
+            if len(self._hamiltonian_filename_Hsys) > 0:
+                mystring += "hamiltonian_file_Hsys= "+str(self._hamiltonian_filename_Hsys)+"\n"
+            if len(self._hamiltonian_filename_Hc) > 0:
+                mystring += "hamiltonian_file_Hc= "+str(self._hamiltonian_filename_Hc)+"\n"
         mystring += "timestepper = "+str(self.timestepper)+ "\n"
         if self.rand_seed is not None and self.rand_seed >= 0:
             mystring += "rand_seed = "+str(int(self.rand_seed))+ "\n"