Quandary
Loading...
Searching...
No Matches
mastereq.hpp
Go to the documentation of this file.
1#include "defs.hpp"
2#include "oscillator.hpp"
3#include "util.hpp"
4#include <optional>
5#include <petscts.h>
6#include <vector>
7#include <assert.h>
8#include <iostream>
9#include "gate.hpp"
10#include "hamiltonianfilereader.hpp"
11
12#pragma once
13
14class MasterEq;
15
16
23typedef struct {
24 MasterEq* mastereq;
25 PetscInt dim;
26 std::vector<size_t> nlevels;
27 IS *isu, *isv;
28 Oscillator** oscil_vec;
29 std::vector<double> crosskerr;
30 std::vector<double> Jkl;
31 std::vector<double> eta;
32 DecoherenceType decoherence_type;
33 bool addT1, addT2;
34 std::vector<double> control_Re;
35 std::vector<double> control_Im;
36 std::vector<Mat> Ac_vec;
37 std::vector<Mat> Bc_vec;
38 Mat *Ad;
39 Mat *Bd;
40 std::vector<Mat> Ad_vec;
41 std::vector<Mat> Bd_vec;
42 std::vector<double> Bd_coeffs; //< Time-dependent coefficients for dipole-dipole coupling matrices: cos(eta_k*t)
43 std::vector<double> Ad_coeffs; //< Time-dependent coefficients for dipole-dipole coupling matrices: sin(eta_k*t)
44 Vec *aux;
45 double time;
46 PetscBool assembled;
47} MatShellCtx;
48
49
58int applyRHS_matfree_1Osc(Mat RHS, Vec x, Vec y);
59int applyRHS_matfree_transpose_1Osc(Mat RHS, Vec x, Vec y);
60int applyRHS_matfree_2Osc(Mat RHS, Vec x, Vec y);
61int applyRHS_matfree_transpose_2Osc(Mat RHS, Vec x, Vec y);
62int applyRHS_matfree_3Osc(Mat RHS, Vec x, Vec y);
63int applyRHS_matfree_transpose_3Osc(Mat RHS, Vec x, Vec y);
64int applyRHS_matfree_4Osc(Mat RHS, Vec x, Vec y);
65int applyRHS_matfree_transpose_4Osc(Mat RHS, Vec x, Vec y);
66int applyRHS_matfree_5Osc(Mat RHS, Vec x, Vec y);
67int applyRHS_matfree_transpose_5Osc(Mat RHS, Vec x, Vec y);
68int applyRHS_sparsemat(Mat RHS, Vec x, Vec y);
69int applyRHS_sparsemat_transpose(Mat RHS, Vec x, Vec y);
70
71
96class MasterEq{
97
98 protected:
99 PetscInt dim;
100 PetscInt dim_rho;
101 PetscInt dim_ess;
102 int noscillators;
103 Oscillator** oscil_vec;
104
105 Mat RHS;
106 MatShellCtx RHSctx;
107
108 std::vector<Mat> Ac_vec; // Vector of constant mats for time-varying control term (real). One for each oscillators.
109 std::vector<Mat> Bc_vec; // Vector of constant mats for time-varying control term (imag). One for each oscillators.
110 Mat Ad, Bd; // Real and imaginary part of constant system matrix
111 std::vector<Mat> Ad_vec; // Vector of constant mats for Dipole-Dipole coupling term in drift Hamiltonian (real)
112 std::vector<Mat> Bd_vec; // Vector of constant mats for Dipole-Dipole coupling term in drift Hamiltonian (imag)
113
114 std::vector<double> crosskerr;
115 std::vector<double> Jkl;
116 std::vector<double> eta;
117 bool addT1, addT2;
118
119 int mpirank_world;
120 int mpirank_petsc;
121 int mpisize_petsc;
122 PetscInt localsize_u;
123 PetscInt ilow;
124 PetscInt iupp;
125
126 IS isu, isv;
127 Vec aux;
128 bool quietmode;
129 std::optional<std::string> hamiltonian_file_Hsys;
130 std::optional<std::string> hamiltonian_file_Hc;
131
132 public:
133 std::vector<size_t> nlevels;
134 std::vector<size_t> nessential;
135 bool usematfree;
136 DecoherenceType decoherence_type;
137
138 public:
139 MasterEq();
140
148 MasterEq(const Config& config, Oscillator** oscil_vec_, bool quietmode=false);
149
150 ~MasterEq();
151
158 void set_RHS_MatMult_operation();
159
163 void initSparseMatSolver();
164
171 Oscillator* getOscillator(const size_t i) { return oscil_vec[i]; }
172
178 size_t getNOscillators() { return noscillators; }
179
185 PetscInt getDim(){ return dim; }
186
192 PetscInt getDimEss(){ return dim_ess; }
193
199 PetscInt getDimRho(){ return dim_rho; }
200
210 int assemble_RHS(const double t);
211
217 Mat getRHS();
218 MatShellCtx* getRHSctx();
219
232 void compute_dRHS_dParams(const double t,const Vec x,const Vec x_bar, const double alpha, Vec grad);
233
239 void setControlAmplitudes(const Vec x);
240
247 double expectedEnergy(const Vec x);
248
255 void population(const Vec x, std::vector<double> &population_com);
256
257 // /* Compute reduced density operator for a sub-system defined by IDs in the oscilIDs vector */
258 // void createReducedDensity(const Vec rho, Vec *reduced, const std::vector<int>& oscilIDs);
259 // /* Derivative of reduced density computation */
260 // void createReducedDensity_diff(Vec rhobar, const Vec reducedbar, const std::vector<int>& oscilIDs);
261};
262
263
283void compute_dRHS_dParams_sparsemat(const double t,const Vec x,const Vec x_bar, const double alpha, Vec grad, std::vector<size_t>& nlevels, IS isu, IS isv, std::vector<Mat>& Ac_vec, std::vector<Mat>& Bc_vec, Vec aux, Oscillator** oscil_vec);
284
301void compute_dRHS_dParams_matfree(const PetscInt dim, const double t,const Vec x,const Vec x_bar, const double alpha, Vec grad, std::vector<size_t>& nlevels, DecoherenceType decoherence_type, Oscillator** oscil_vec);
302
303
304
305// Inline functions for the Matrix-free RHS application
306
314inline double H_detune(const double detuning0, const int a) {
315 return detuning0*a;
316};
317
325inline double H_selfkerr(const double xi0, const int a) {
326 return - xi0 / 2.0 * a * (a-1);
327};
328
337inline double L2(const double dephase0, const int i0, const int i0p){
338 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) );
339};
340
349inline double L1diag(const double decay0, const int i0, const int i0p){
350 return - decay0 / 2.0 * ( i0 + i0p ) ;
351};
352
361inline int TensorGetIndex(const int nlevels0, const int i0, const int i0p){
362 return i0 + nlevels0 * i0p;
363};
364
374inline double H_detune(const double detuning0, const double detuning1, const int a, const int b) {
375 return detuning0*a + detuning1*b;
376};
377
387inline double H_selfkerr(const double xi0, const double xi1, const int a, const int b) {
388 return - xi0 / 2.0 * a * (a-1) - xi1 / 2.0 * b * (b-1);
389};
390
399inline double H_crosskerr(const double xi01, const int a, const int b) {
400 return - xi01 * a * b;
401};
402
414inline double L2(const double dephase0, const double dephase1, const int i0, const int i1, const int i0p, const int i1p){
415 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) ) + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) );
416};
417
429inline double L1diag(const double decay0, const double decay1, const int i0, const int i1, const int i0p, const int i1p){
430 return - decay0 / 2.0 * ( i0 + i0p ) - decay1 / 2.0 * ( i1 + i1p );
431};
432
444inline int TensorGetIndex(const int nlevels0, const int nlevels1,const int i0, const int i1, const int i0p, const int i1p){
445 return i0*nlevels1 + i1 + (nlevels0 * nlevels1) * ( i0p * nlevels1 + i1p);
446};
447
448
449// Inlines for Matrix-free RHS for 3 oscillators
450// See documentation for 1-2 oscillator functions above for parameter descriptions
451inline double H_detune(const double detuning0, const double detuning1, const double detuning2, const int i0, const int i1, const int i2) {
452 return detuning0*i0 + detuning1*i1 + detuning2*i2;
453};
454inline double H_selfkerr(const double xi0, const double xi1, const double xi2, const int i0, const int i1, const int i2) {
455 return - xi0 / 2.0 * i0 * (i0-1) - xi1 / 2.0 * i1 * (i1-1) - xi2 / 2.0 * i2 * (i2-1);
456};
457inline double H_crosskerr(const double xi01, const double xi02, const double xi12, const int i0, const int i1, const int i2) {
458 return - xi01 * i0 * i1 - xi02 * i0 * i2 - xi12 * i1 * i2;
459};
460inline double L2(const double dephase0, const double dephase1, const double dephase2, const int i0, const int i1, const int i2, const int i0p, const int i1p, const int i2p){
461 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) )
462 + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) )
463 + dephase2 * ( i2*i2p - 1./2. * (i2*i2 + i2p*i2p) );
464};
465inline double L1diag(const double decay0, const double decay1, const double decay2, const int i0, const int i1, const int i2, const int i0p, const int i1p, const int i2p){
466 return - decay0 / 2.0 * ( i0 + i0p )
467 - decay1 / 2.0 * ( i1 + i1p )
468 - decay2 / 2.0 * ( i2 + i2p );
469};
470inline int TensorGetIndex(const int nlevels0, const int nlevels1, const int nlevels2, const int i0, const int i1, const int i2, const int i0p, const int i1p, const int i2p){
471 return i0*nlevels1*nlevels2 + i1*nlevels2 + i2 + (nlevels0 * nlevels1 * nlevels2) * ( i0p * nlevels1*nlevels2 + i1p*nlevels2 + i2p);
472};
473
474
475// Matrix-free solver inlines for 4 oscillators
476// See documentation for 1-2 oscillator functions above for parameter descriptions
477inline double H_detune(const double detuning0, const double detuning1, const double detuning2, const double detuning3, const int i0, const int i1, const int i2, const int i3) {
478 return detuning0*i0 + detuning1*i1 + detuning2*i2 + detuning3*i3;
479};
480inline double H_selfkerr(const double xi0, const double xi1, const double xi2, const double xi3, const int i0, const int i1, const int i2, const int i3) {
481 return - xi0 / 2.0 * i0 * (i0-1) - xi1 / 2.0 * i1 * (i1-1) - xi2 / 2.0 * i2 * (i2-1) - xi3/2.0 * i3 * (i3-1);
482};
483inline double H_crosskerr(const double xi01, const double xi02, const double xi03, const double xi12, const double xi13, const double xi23, const int i0, const int i1, const int i2, const int i3) {
484 return - xi01 * i0 * i1 - xi02 * i0 * i2 - xi03*i0*i3 - xi12 * i1 * i2 - xi13*i1*i3 - xi23*i2*i3;
485};
486inline double L2(const double dephase0, const double dephase1, const double dephase2, const double dephase3, const int i0, const int i1, const int i2, const int i3, const int i0p, const int i1p, const int i2p, const int i3p){
487 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) )
488 + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) )
489 + dephase2 * ( i2*i2p - 1./2. * (i2*i2 + i2p*i2p) )
490 + dephase3 * ( i3*i3p - 1./2. * (i3*i3 + i3p*i3p) );
491};
492inline double L1diag(const double decay0, const double decay1, const double decay2, const double decay3, const int i0, const int i1, const int i2, const int i3, const int i0p, const int i1p, const int i2p, const int i3p){
493 return - decay0 / 2.0 * ( i0 + i0p )
494 - decay1 / 2.0 * ( i1 + i1p )
495 - decay2 / 2.0 * ( i2 + i2p )
496 - decay3 / 2.0 * ( i3 + i3p );
497};
498inline int TensorGetIndex(const int nlevels0, const int nlevels1, const int nlevels2, const int nlevels3, const int i0, const int i1, const int i2, const int i3, const int i0p, const int i1p, const int i2p, const int i3p){
499 return i0*nlevels1*nlevels2*nlevels3 + i1*nlevels2*nlevels3 + i2*nlevels3 + i3 + (nlevels0 * nlevels1 * nlevels2 * nlevels3) * ( i0p * nlevels1*nlevels2*nlevels3 + i1p*nlevels2*nlevels3 + i2p*nlevels3 + i3p);
500}
501
502// Matrix-free solver inlines for 5 oscillators
503// See documentation for 1-2 oscillator functions above for parameter descriptions
504inline double H_detune(const double detuning0, const double detuning1, const double detuning2, const double detuning3, const double detuning4, const int i0, const int i1, const int i2, const int i3, const int i4) {
505 return detuning0*i0 + detuning1*i1 + detuning2*i2 + detuning3*i3 + detuning4*i4;
506};
507inline double H_selfkerr(const double xi0, const double xi1, const double xi2, const double xi3, const double xi4, const int i0, const int i1, const int i2, const int i3, const int i4) {
508 return - xi0 / 2.0 * i0 * (i0-1) - xi1 / 2.0 * i1 * (i1-1) - xi2 / 2.0 * i2 * (i2-1) - xi3/2.0 * i3 * (i3-1) - xi4/2.0 * i4 * (i4-1);
509};
510inline double H_crosskerr(const double xi01, const double xi02, const double xi03, const double xi04, const double xi12, const double xi13, const double xi14, const double xi23, const double xi24, const double xi34, const int i0, const int i1, const int i2, const int i3, const int i4) {
511 return - xi01 * i0 * i1 - xi02 * i0 * i2 - xi03*i0*i3 - xi04*i0*i4 - xi12 * i1 * i2 - xi13*i1*i3 - xi14*i1*i4 - xi23*i2*i3 - xi24*i2*i4 - xi34*i3*i4;
512};
513inline double L2(const double dephase0, const double dephase1, const double dephase2, const double dephase3, const double dephase4, const int i0, const int i1, const int i2, const int i3, const int i4, const int i0p, const int i1p, const int i2p, const int i3p, const int i4p){
514 return dephase0 * ( i0*i0p - 1./2. * (i0*i0 + i0p*i0p) )
515 + dephase1 * ( i1*i1p - 1./2. * (i1*i1 + i1p*i1p) )
516 + dephase2 * ( i2*i2p - 1./2. * (i2*i2 + i2p*i2p) )
517 + dephase3 * ( i3*i3p - 1./2. * (i3*i3 + i3p*i3p) )
518 + dephase4 * ( i4*i4p - 1./2. * (i4*i4 + i4p*i4p) );
519};
520inline double L1diag(const double decay0, const double decay1, const double decay2, const double decay3, const double decay4, const int i0, const int i1, const int i2, const int i3, const int i4, const int i0p, const int i1p, const int i2p, const int i3p, const int i4p){
521 return - decay0 / 2.0 * ( i0 + i0p )
522 - decay1 / 2.0 * ( i1 + i1p )
523 - decay2 / 2.0 * ( i2 + i2p )
524 - decay3 / 2.0 * ( i3 + i3p )
525 - decay4 / 2.0 * ( i4 + i4p );
526};
527inline int TensorGetIndex(const int nlevels0, const int nlevels1, const int nlevels2, const int nlevels3, const int nlevels4, const int i0, const int i1, const int i2, const int i3, const int i4, const int i0p, const int i1p, const int i2p, const int i3p, const int i4p){
528 return i0*nlevels1*nlevels2*nlevels3*nlevels4 + i1*nlevels2*nlevels3*nlevels4 + i2*nlevels3*nlevels4 + i3*nlevels4 + i4 + (nlevels0 * nlevels1 * nlevels2 * nlevels3*nlevels4) * ( i0p * nlevels1*nlevels2*nlevels3*nlevels4 + i1p*nlevels2*nlevels3*nlevels4 + i2p*nlevels3*nlevels4+ i3p*nlevels4 + i4p);
529}
530
531
551inline void dRHSdp_getcoeffs(const PetscInt dim, const int it, const int n, const int np, const int i, const int ip, const int stridei, const int strideip, const double* xptr, double* res_p_re, double* res_p_im, double* res_q_re, double* res_q_im) {
552
553 *res_p_re = 0.0;
554 *res_p_im = 0.0;
555 *res_q_re = 0.0;
556 *res_q_im = 0.0;
557
558 /* ik+1..,ik'.. term */
559 if (i < n-1) {
560 int itx = it + stridei;
561 double xre = xptr[itx];
562 double xim = xptr[itx+dim];
563 double sq = sqrt(i + 1);
564 *res_p_re += sq * xim;
565 *res_p_im += - sq * xre;
566 *res_q_re += sq * xre;
567 *res_q_im += sq * xim;
568 }
569 /* \rho(ik..,ik'+1..) */
570 if (ip < np-1) {
571 int itx = it + strideip;
572 double xre = xptr[itx];
573 double xim = xptr[itx+dim];
574 double sq = sqrt(ip + 1);
575 *res_p_re += - sq * xim;
576 *res_p_im += + sq * xre;
577 *res_q_re += sq * xre;
578 *res_q_im += sq * xim;
579 }
580 /* \rho(ik-1..,ik'..) */
581 if (i > 0) {
582 int itx = it - stridei;
583 double xre = xptr[itx];
584 double xim = xptr[itx+dim];
585 double sq = sqrt(i);
586 *res_p_re += + sq * xim;
587 *res_p_im += - sq * xre;
588 *res_q_re += - sq * xre;
589 *res_q_im += - sq * xim;
590 }
591 /* \rho(ik..,ik'-1..) */
592 if (ip > 0) {
593 int itx = it - strideip;
594 double xre = xptr[itx];
595 double xim = xptr[itx+dim];
596 double sq = sqrt(ip);
597 *res_p_re += - sq * xim;
598 *res_p_im += + sq * xre;
599 *res_q_re += - sq * xre;
600 *res_q_im += - sq * xim;
601 }
602}
603
630inline void Jkl_coupling(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double* xptr, const double Jij, const double cosij, const double sinij, double* yre, double* yim) {
631 if (fabs(Jij)>1e-10) {
632 // 1) J_kl (-icos + sin) * ρ_{E−k+l i, i′}
633 if (i > 0 && j < nj-1) {
634 int itx = it - stridei + stridej;
635 double xre = xptr[itx];
636 double xim = xptr[itx+dim];
637 double sq = sqrt(i * (j + 1));
638 // sin u + cos v + i ( -cos u + sin v)
639 *yre += Jij * sq * ( cosij * xim + sinij * xre);
640 *yim += Jij * sq * ( - cosij * xre + sinij * xim);
641 }
642 // 2) J_kl (−icos − sin)sqrt(il*(ik +1)) ρ_{E+k−li,i′}
643 if (i < ni-1 && j > 0) {
644 int itx = it + stridei - stridej; // E+k-l i, i'
645 double xre = xptr[itx];
646 double xim = xptr[itx+dim];
647 double sq = sqrt(j * (i + 1)); // sqrt( il*(ik+1))
648 // -sin u + cos v + i (-cos u - sin v)
649 *yre += Jij * sq * ( cosij * xim - sinij * xre);
650 *yim += Jij * sq * ( - cosij * xre - sinij * xim);
651 }
652 // 3) J_kl ( icos + sin)sqrt(ik'*(il' +1)) ρ_{i,E-k+li'}
653 if (ip > 0 && jp < njp-1) {
654 int itx = it - strideip + stridejp; // i, E-k+l i'
655 double xre = xptr[itx];
656 double xim = xptr[itx+dim];
657 double sq = sqrt(ip * (jp + 1)); // sqrt( ik'*(il'+1))
658 // sin u - cos v + i ( cos u + sin v)
659 *yre += Jij * sq * ( - cosij * xim + sinij * xre);
660 *yim += Jij * sq * ( cosij * xre + sinij * xim);
661 }
662 // 4) J_kl ( icos - sin)sqrt(il'*(ik' +1)) ρ_{i,E+k-li'}
663 if (ip < nip-1 && jp > 0) {
664 int itx = it + strideip - stridejp; // i, E+k-l i'
665 double xre = xptr[itx];
666 double xim = xptr[itx+dim];
667 double sq = sqrt(jp * (ip + 1)); // sqrt( il'*(ik'+1))
668 // - sin u - cos v + i ( cos u - sin v)
669 *yre += Jij * sq * ( - cosij * xim - sinij * xre);
670 *yim += Jij * sq * ( cosij * xre - sinij * xim);
671 }
672 }
673}
674
699inline void Jkl_coupling_T(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double* xptr, const double Jij, const double cosij, const double sinij, double* yre, double* yim) {
700 if (fabs(Jij)>1e-10) {
701 // 1) [...] * \bar y_{E+k-l i, i′}
702 if (i < ni-1 && j > 0) {
703 int itx = it + stridei - stridej;
704 double xre = xptr[itx];
705 double xim = xptr[itx+dim];
706 double sq = sqrt(j * (i + 1));
707 *yre += Jij * sq * ( - cosij * xim + sinij * xre);
708 *yim += Jij * sq * ( + cosij * xre + sinij * xim);
709 }
710 // 2) J_kl (−icos − sin)sqrt(ik*(il +1)) \bar y_{E-k+li,i′}
711 if (i > 0 && j < nj-1) {
712 int itx = it - stridei + stridej; // E-k+l i, i'
713 double xre = xptr[itx];
714 double xim = xptr[itx+dim];
715 double sq = sqrt(i * (j + 1)); // sqrt( ik*(il+1))
716 *yre += Jij * sq * ( - cosij * xim - sinij * xre);
717 *yim += Jij * sq * ( + cosij * xre - sinij * xim);
718 }
719 // 3) J_kl ( icos + sin)sqrt(il'*(ik' +1)) \bar y_{i,E+k-li'}
720 if (ip < nip-1 && jp > 0) {
721 int itx = it + strideip - stridejp; // i, E+k-l i'
722 double xre = xptr[itx];
723 double xim = xptr[itx+dim];
724 double sq = sqrt(jp * (ip + 1)); // sqrt( il'*(ik'+1))
725 *yre += Jij * sq * ( cosij * xim + sinij * xre);
726 *yim += Jij * sq * ( - cosij * xre + sinij * xim);
727 }
728 // 4) J_kl ( icos - sin)sqrt(ik'*(il' +1)) \bar y_{i,E-k+li'}
729 if (ip > 0 && jp < njp-1) {
730 int itx = it - strideip + stridejp; // i, E-k+l i'
731 double xre = xptr[itx];
732 double xim = xptr[itx+dim];
733 double sq = sqrt(ip * (jp + 1)); // sqrt( ik'*(il'+1))
734 *yre += Jij * sq * ( cosij * xim - sinij * xre);
735 *yim += Jij * sq * ( - cosij * xre - sinij * xim);
736 }
737 }
738
739}
740
756inline void L1decay(const PetscInt dim, const int it, const int n, const int i, const int ip, const int stridei, const int strideip, const double* xptr, const double decayi, double* yre, double* yim){
757 if (fabs(decayi) > 1e-12) {
758 if (i < n-1 && ip < n-1) {
759 double l1off = decayi * sqrt((i+1)*(ip+1));
760 int itx = it + stridei + strideip;
761 double xre = xptr[itx];
762 double xim = xptr[itx+dim];
763 *yre += l1off * xre;
764 *yim += l1off * xim;
765 }
766 }
767}
768
769
784inline void L1decay_T(const PetscInt dim, const int it, const int i, const int ip, const int stridei, const int strideip, const double* xptr, const double decayi, double* yre, double* yim){
785 if (fabs(decayi) > 1e-12) {
786 if (i > 0 && ip > 0) {
787 double l1off = decayi * sqrt(i*ip);
788 int itx = it - stridei - strideip;
789 double xre = xptr[itx];
790 double xim = xptr[itx+dim];
791 *yre += l1off * xre;
792 *yim += l1off * xim;
793 }
794 }
795}
796
816inline void control(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double* xptr, const double pt, const double qt, double* yre, double* yim){
817 /* \rho(ik+1..,ik'..) term */
818 if (i < n-1) {
819 int itx = it + stridei;
820 double xre = xptr[itx];
821 double xim = xptr[itx+dim];
822 double sq = sqrt(i + 1);
823 *yre += sq * ( pt * xim + qt * xre);
824 *yim += sq * ( - pt * xre + qt * xim);
825 }
826 /* \rho(ik..,ik'+1..) */
827 if (ip < np-1) {
828 int itx = it + strideip;
829 double xre = xptr[itx];
830 double xim = xptr[itx+dim];
831 double sq = sqrt(ip + 1);
832 *yre += sq * ( -pt * xim + qt * xre);
833 *yim += sq * ( pt * xre + qt * xim);
834 }
835 /* \rho(ik-1..,ik'..) */
836 if (i > 0) {
837 int itx = it - stridei;
838 double xre = xptr[itx];
839 double xim = xptr[itx+dim];
840 double sq = sqrt(i);
841 *yre += sq * ( pt * xim - qt * xre);
842 *yim += sq * (- pt * xre - qt * xim);
843 }
844 /* \rho(ik..,ik'-1..) */
845 if (ip > 0) {
846 int itx = it - strideip;
847 double xre = xptr[itx];
848 double xim = xptr[itx+dim];
849 double sq = sqrt(ip);
850 *yre += sq * (- pt * xim - qt * xre);
851 *yim += sq * ( pt * xre - qt * xim);
852 }
853}
854
855
873inline void control_T(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double* xptr, const double pt, const double qt, double* yre, double* yim){
874 /* \rho(ik+1..,ik'..) term */
875 if (i > 0) {
876 int itx = it - stridei;
877 double xre = xptr[itx];
878 double xim = xptr[itx+dim];
879 double sq = sqrt(i);
880 *yre += sq * ( - pt * xim + qt * xre);
881 *yim += sq * ( pt * xre + qt * xim);
882 }
883 /* \rho(ik..,ik'+1..) */
884 if (ip > 0) {
885 int itx = it - strideip;
886 double xre = xptr[itx];
887 double xim = xptr[itx+dim];
888 double sq = sqrt(ip);
889 *yre += sq * ( pt * xim + qt * xre);
890 *yim += sq * ( -pt * xre + qt * xim);
891 }
892 /* \rho(ik-1..,ik'..) */
893 if (i < n-1) {
894 int itx = it + stridei;
895 double xre = xptr[itx];
896 double xim = xptr[itx+dim];
897 double sq = sqrt(i+1);
898 *yre += sq * (- pt * xim - qt * xre);
899 *yim += sq * ( pt * xre - qt * xim);
900 }
901 /* \rho(ik..,ik'-1..) */
902 if (ip < np-1) {
903 int itx = it + strideip;
904 double xre = xptr[itx];
905 double xim = xptr[itx+dim];
906 double sq = sqrt(ip+1);
907 *yre += sq * (+ pt * xim - qt * xre);
908 *yim += sq * (- pt * xre - qt * xim);
909 }
910}
Config
Configuration class containing all validated settings.
Definition config.hpp:56
MasterEq
Implementation of the real-valued right-hand-side (RHS) system matrix of the quantum dynamical equati...
Definition mastereq.hpp:96
MasterEq::addT1
bool addT1
Definition mastereq.hpp:117
MasterEq::addT2
bool addT2
Flags for including T1 decay and T2 dephasing Lindblad operators.
Definition mastereq.hpp:117
MasterEq::Ac_vec
std::vector< Mat > Ac_vec
Definition mastereq.hpp:108
MasterEq::decoherence_type
DecoherenceType decoherence_type
Type of Lindblad operators to include (NONE means Schroedinger equation)
Definition mastereq.hpp:136
MasterEq::getRHS
Mat getRHS()
Retrieves the right-hand-side system matrix.
Definition mastereq.cpp:703
MasterEq::nlevels
std::vector< size_t > nlevels
Number of levels per oscillator.
Definition mastereq.hpp:133
MasterEq::mpirank_petsc
int mpirank_petsc
Rank of PETSc's communicator.
Definition mastereq.hpp:120
MasterEq::MasterEq
MasterEq()
Definition mastereq.cpp:3
MasterEq::Bd
Mat Bd
Definition mastereq.hpp:110
MasterEq::usematfree
bool usematfree
Flag for using matrix-free solver.
Definition mastereq.hpp:135
MasterEq::dim_rho
PetscInt dim_rho
Dimension of Hilbert space = N.
Definition mastereq.hpp:100
MasterEq::quietmode
bool quietmode
Flag for quiet mode operation.
Definition mastereq.hpp:128
MasterEq::nessential
std::vector< size_t > nessential
Number of essential levels per oscillator.
Definition mastereq.hpp:134
MasterEq::crosskerr
std::vector< double > crosskerr
Cross-Kerr coefficients (rad/time) for ZZ-coupling .
Definition mastereq.hpp:114
MasterEq::isu
IS isu
Definition mastereq.hpp:126
MasterEq::hamiltonian_file_Hc
std::optional< std::string > hamiltonian_file_Hc
Filename if a custom control Hamiltonians are read from file.
Definition mastereq.hpp:130
MasterEq::isv
IS isv
Vector strides for accessing real and imaginary parts u=Re(x), v=Im(x)
Definition mastereq.hpp:126
MasterEq::localsize_u
PetscInt localsize_u
Size of local sub vector u or v in state x=[u,v].
Definition mastereq.hpp:122
MasterEq::getDim
PetscInt getDim()
Retrieves the dimension of the vectorized system.
Definition mastereq.hpp:185
MasterEq::getNOscillators
size_t getNOscillators()
Retrieves the number of oscillators in the system.
Definition mastereq.hpp:178
MasterEq::getDimRho
PetscInt getDimRho()
Retrieves the dimension of the density matrix.
Definition mastereq.hpp:199
MasterEq::RHSctx
MatShellCtx RHSctx
Context containing data for applying the RHS to a state.
Definition mastereq.hpp:106
MasterEq::oscil_vec
Oscillator ** oscil_vec
Array of pointers to oscillator objects.
Definition mastereq.hpp:103
MasterEq::Bd_vec
std::vector< Mat > Bd_vec
Definition mastereq.hpp:112
MasterEq::hamiltonian_file_Hsys
std::optional< std::string > hamiltonian_file_Hsys
Filename if a custom system Hamiltonian is read from file.
Definition mastereq.hpp:129
MasterEq::Ad
Mat Ad
Definition mastereq.hpp:110
MasterEq::getOscillator
Oscillator * getOscillator(const size_t i)
Retrieves the i-th oscillator.
Definition mastereq.hpp:171
MasterEq::eta
std::vector< double > eta
Frequency differences in rotating frame (rad/time) for dipole-dipole coupling.
Definition mastereq.hpp:116
MasterEq::population
void population(const Vec x, std::vector< double > &population_com)
Computes population of the full composite system.
Definition mastereq.cpp:2963
MasterEq::set_RHS_MatMult_operation
void set_RHS_MatMult_operation()
Pass MatMult operations for applying the RHS to PETSc.
Definition mastereq.cpp:737
MasterEq::getRHSctx
MatShellCtx * getRHSctx()
Definition mastereq.cpp:704
MasterEq::ilow
PetscInt ilow
First index of the local sub vector u,v.
Definition mastereq.hpp:123
MasterEq::noscillators
int noscillators
Number of oscillators in the system.
Definition mastereq.hpp:102
MasterEq::initSparseMatSolver
void initSparseMatSolver()
Initializes matrices needed for the sparse matrix solver.
Definition mastereq.cpp:214
MasterEq::assemble_RHS
int assemble_RHS(const double t)
Assembles the real-valued system matrix (RHS) at time t.
Definition mastereq.cpp:679
MasterEq::expectedEnergy
double expectedEnergy(const Vec x)
Computes expected energy of the full composite system.
Definition mastereq.cpp:2924
MasterEq::Ad_vec
std::vector< Mat > Ad_vec
Definition mastereq.hpp:111
MasterEq::RHS
Mat RHS
MatShell for real-valued, vectorized system matrix (size 2N^2 x 2N^2)
Definition mastereq.hpp:105
MasterEq::mpisize_petsc
int mpisize_petsc
Size of PETSc's communicator.
Definition mastereq.hpp:121
MasterEq::Jkl
std::vector< double > Jkl
Dipole-dipole coupling coefficients (rad/time), multiplies .
Definition mastereq.hpp:115
MasterEq::mpirank_world
int mpirank_world
Rank of global MPI communicator.
Definition mastereq.hpp:119
MasterEq::aux
Vec aux
Auxiliary vector for computations.
Definition mastereq.hpp:127
MasterEq::dim
PetscInt dim
Dimension of full vectorized system: N^2 if Lindblad, N if Schroedinger.
Definition mastereq.hpp:99
MasterEq::getDimEss
PetscInt getDimEss()
Retrieves the dimension of the essential level system.
Definition mastereq.hpp:192
MasterEq::dim_ess
PetscInt dim_ess
Dimension of essential levels = N_e.
Definition mastereq.hpp:101
MasterEq::iupp
PetscInt iupp
Last index (+1) of the local sub vector u,v.
Definition mastereq.hpp:124
MasterEq::Bc_vec
std::vector< Mat > Bc_vec
Definition mastereq.hpp:109
MasterEq::setControlAmplitudes
void setControlAmplitudes(const Vec x)
Pass control parameters from global design vector to each oscillator.
Definition mastereq.cpp:717
MasterEq::~MasterEq
~MasterEq()
Definition mastereq.cpp:180
MasterEq::compute_dRHS_dParams
void compute_dRHS_dParams(const double t, const Vec x, const Vec x_bar, const double alpha, Vec grad)
Computes gradient of RHS with respect to control parameters.
Definition mastereq.cpp:707
Oscillator
Quantum oscillator (multi-level qubit) with control capabilities.
Definition oscillator.hpp:37
defs.hpp
Core type definitions and enumerations for Quandary quantum optimal control.
DecoherenceType
DecoherenceType
Available Lindblad operator types for open quantum systems, or NONE for closed quantum systems.
Definition defs.hpp:31
compute_dRHS_dParams_sparsemat
void compute_dRHS_dParams_sparsemat(const double t, const Vec x, const Vec x_bar, const double alpha, Vec grad, std::vector< size_t > &nlevels, IS isu, IS isv, std::vector< Mat > &Ac_vec, std::vector< Mat > &Bc_vec, Vec aux, Oscillator **oscil_vec)
: Sparse-matrix version to compute gradient of RHS with respect to parameters
Definition mastereq.cpp:952
Jkl_coupling_T
void Jkl_coupling_T(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double *xptr, const double Jij, const double cosij, const double sinij, double *yre, double *yim)
Transpose of dipole-dipole coupling for adjoint computations.
Definition mastereq.hpp:699
applyRHS_sparsemat
int applyRHS_sparsemat(Mat RHS, Vec x, Vec y)
Sparse matrix MatMult.
Definition mastereq.cpp:770
compute_dRHS_dParams_matfree
void compute_dRHS_dParams_matfree(const PetscInt dim, const double t, const Vec x, const Vec x_bar, const double alpha, Vec grad, std::vector< size_t > &nlevels, DecoherenceType decoherence_type, Oscillator **oscil_vec)
: Matrix free version to compute gradient of RHS with respect to parameters
Definition mastereq.cpp:997
applyRHS_matfree_1Osc
int applyRHS_matfree_1Osc(Mat RHS, Vec x, Vec y)
Matrix-vector products to apply the RHS system matrix to a state vector.
Definition mastereq.cpp:3004
L2
double L2(const double dephase0, const int i0, const int i0p)
Inline for Matrix-free RHS to compute L2 dephasing coefficient for 1 oscillator.
Definition mastereq.hpp:337
applyRHS_matfree_transpose_1Osc
int applyRHS_matfree_transpose_1Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 1 oscillator.
Definition mastereq.cpp:3033
applyRHS_matfree_transpose_5Osc
int applyRHS_matfree_transpose_5Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 5 oscillators.
Definition mastereq.cpp:3245
applyRHS_matfree_transpose_2Osc
int applyRHS_matfree_transpose_2Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 2 oscillators.
Definition mastereq.cpp:3090
applyRHS_matfree_4Osc
int applyRHS_matfree_4Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 4 oscillators.
Definition mastereq.cpp:3168
applyRHS_matfree_3Osc
int applyRHS_matfree_3Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 3 oscillators.
Definition mastereq.cpp:3118
dRHSdp_getcoeffs
void dRHSdp_getcoeffs(const PetscInt dim, const int it, const int n, const int np, const int i, const int ip, const int stridei, const int strideip, const double *xptr, double *res_p_re, double *res_p_im, double *res_q_re, double *res_q_im)
Inline for Matrix-free RHS for gradient updates.
Definition mastereq.hpp:551
Jkl_coupling
void Jkl_coupling(const PetscInt dim, const int it, const int ni, const int nj, const int nip, const int njp, const int i, const int ip, const int j, const int jp, const int stridei, const int strideip, const int stridej, const int stridejp, const double *xptr, const double Jij, const double cosij, const double sinij, double *yre, double *yim)
Matrix-free solver inline for dipole-dipole coupling between oscillators.
Definition mastereq.hpp:630
L1decay
void L1decay(const PetscInt dim, const int it, const int n, const int i, const int ip, const int stridei, const int strideip, const double *xptr, const double decayi, double *yre, double *yim)
Matrix-free solver inline for off-diagonal L1 decay term.
Definition mastereq.hpp:756
H_crosskerr
double H_crosskerr(const double xi01, const int a, const int b)
Inline for Matrix-free RHS, Computes cross-Kerr coupling between 2 oscillators.
Definition mastereq.hpp:399
TensorGetIndex
int TensorGetIndex(const int nlevels0, const int i0, const int i0p)
Inline for Matrix-free RHS to compute tensor product index for 1 oscillator system.
Definition mastereq.hpp:361
control
void control(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double *xptr, const double pt, const double qt, double *yre, double *yim)
Matrix-free solver inline for control terms.
Definition mastereq.hpp:816
applyRHS_matfree_5Osc
int applyRHS_matfree_5Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 5 oscillators.
Definition mastereq.cpp:3219
H_selfkerr
double H_selfkerr(const double xi0, const int a)
Inline for Matrix-free RHS to Compute self-Kerr coefficient for 1 oscillator.
Definition mastereq.hpp:325
applyRHS_matfree_2Osc
int applyRHS_matfree_2Osc(Mat RHS, Vec x, Vec y)
Matrix-free MatMult for 2 oscillators.
Definition mastereq.cpp:3060
applyRHS_matfree_transpose_4Osc
int applyRHS_matfree_transpose_4Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 4 oscillators.
Definition mastereq.cpp:3194
control_T
void control_T(const PetscInt dim, const int it, const int n, const int i, const int np, const int ip, const int stridei, const int strideip, const double *xptr, const double pt, const double qt, double *yre, double *yim)
Matrix-free Transpose of control terms for adjoint computations.
Definition mastereq.hpp:873
L1diag
double L1diag(const double decay0, const int i0, const int i0p)
Inline for Matrix-free RHS to compute L1 decay coefficient for 1 oscillator.
Definition mastereq.hpp:349
applyRHS_matfree_transpose_3Osc
int applyRHS_matfree_transpose_3Osc(Mat RHS, Vec x, Vec y)
Transpose matrix-free MatMult for 3 oscillators.
Definition mastereq.cpp:3144
applyRHS_sparsemat_transpose
int applyRHS_sparsemat_transpose(Mat RHS, Vec x, Vec y)
Transpose sparse matrix MatMult.
Definition mastereq.cpp:861
H_detune
double H_detune(const double detuning0, const int a)
Inline for Matrix-free RHS to Compute detuning for 1 oscillator.
Definition mastereq.hpp:314
L1decay_T
void L1decay_T(const PetscInt dim, const int it, const int i, const int ip, const int stridei, const int strideip, const double *xptr, const double decayi, double *yre, double *yim)
Matrix-free inline Transpose of off-diagonal L1 decay for adjoint computations.
Definition mastereq.hpp:784
MatShellCtx
Matrix shell context containing data needed for applying the right-hand-side (RHS) system matrix to a...
Definition mastereq.hpp:23
MatShellCtx::Ad_vec
std::vector< Mat > Ad_vec
Vector of real parts of dipole-dipole coupling system matrices.
Definition mastereq.hpp:40
MatShellCtx::mastereq
MasterEq * mastereq
Owning MasterEq instance for callback dispatch.
Definition mastereq.hpp:24
MatShellCtx::oscil_vec
Oscillator ** oscil_vec
Array of pointers to the oscillators.
Definition mastereq.hpp:28
MatShellCtx::aux
Vec * aux
Auxiliary vector for computations.
Definition mastereq.hpp:44
MatShellCtx::Jkl
std::vector< double > Jkl
Dipole-dipole coupling strength.
Definition mastereq.hpp:30
MatShellCtx::eta
std::vector< double > eta
Frequency differences of the rotating frames.
Definition mastereq.hpp:31
MatShellCtx::nlevels
std::vector< size_t > nlevels
Number of levels per oscillator.
Definition mastereq.hpp:26
MatShellCtx::crosskerr
std::vector< double > crosskerr
Cross-Kerr coupling coefficients.
Definition mastereq.hpp:29
MatShellCtx::time
double time
Current time.
Definition mastereq.hpp:45
MatShellCtx::Ad
Mat * Ad
Real parts of time-independent system matrix.
Definition mastereq.hpp:38
MatShellCtx::control_Im
std::vector< double > control_Im
Imaginary parts of control pulse .
Definition mastereq.hpp:35
MatShellCtx::dim
PetscInt dim
Dimension of full vectorized system: N^2 if Lindblad, N if Schroedinger.
Definition mastereq.hpp:25
MatShellCtx::control_Re
std::vector< double > control_Re
Real parts of control pulse .
Definition mastereq.hpp:34
MatShellCtx::Bd_coeffs
std::vector< double > Bd_coeffs
Definition mastereq.hpp:42
MatShellCtx::Ac_vec
std::vector< Mat > Ac_vec
Vector of real parts of control matrices per oscillator.
Definition mastereq.hpp:36
MatShellCtx::isu
IS * isu
Definition mastereq.hpp:27
MatShellCtx::Bd
Mat * Bd
Imaginary parts of time-independent system matrix.
Definition mastereq.hpp:39
MatShellCtx::Ad_coeffs
std::vector< double > Ad_coeffs
Definition mastereq.hpp:43
MatShellCtx::decoherence_type
DecoherenceType decoherence_type
Type of Lindblad operators to include.
Definition mastereq.hpp:32
MatShellCtx::addT1
bool addT1
Definition mastereq.hpp:33
MatShellCtx::assembled
PetscBool assembled
Definition mastereq.hpp:46
MatShellCtx::Bd_vec
std::vector< Mat > Bd_vec
Vector of imaginary parts of dipole-dipole coupling system matrices.
Definition mastereq.hpp:41
MatShellCtx::Bc_vec
std::vector< Mat > Bc_vec
Vector of imaginary parts of control matrices per oscillator.
Definition mastereq.hpp:37
Generated by doxygen 1.9.8