Hi, interesting set of questions but it's hard to give a totally general answer, as it could depend on many things.
The first and most important question is: have you done careful profiling/timing of the code to be sure about which part is taking the most time? Sometimes it can be surprising as to which part is taking the longest if you measure.
Assuming you've now done that, then for the two likely slowest parts here are some of the things you can try:
(0) you can just try lowering the accuracy you are asking for (larger cutoff, limiting maxm) or using a bigger time step, and see if you still get acceptable results
(1) if a slow part is the decomposition of the local piece of the wavefunction after changing it (the call to svdBond), then you can make extra sure that the method used to factor the wavefunction tensor is the "denmatDecomp" method and not the "svd" method. The .svdBond method (confusingly) switches between these depending on the particular accuracy parameters you set; if you are setting a rather small cutoff, the code may be using the more accurate svd routine versus the faster denmatDecomp routine.
(2) if a slow part is exponentiating your operators, and you can justify separately exponentiating your diagonal operators without causing too much of an additional Trotter / time-step error (there may not be any extra error; just something to make sure to think about) then yes there would be an analytical expression for the exponential of a diagonal operator of course, so you can just directly set the elements of the associated ITensor/IQTensor you want to construct by using the .set method.
Regarding the .set method, you may be interested to know we just added a more convenient way to set multiple elements:
auto T = ITensor(i,j,k);
T.order(k,j,i); //ensures that the order of the indices is k,j,i
T.set(2,1,3,0.5235); //set the k=2,j=1,i=3 element to 0.5235