|
6 | 6 | %% FEYN DYN, VERSION 2013.11.28 |
7 | 7 | %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
8 | 8 |
|
9 | | -function [rho,elapsedTime]=FeynDyn(finalPoint,deltaKmax,totalT,rho,H,systemCouplingMatrix,w,dw,J,temperature,wholeDensityMatrixOrJustDiagonals,allPointsORjustFinalPoint,cpuORgpu) |
| 9 | +function [rho,elapsedTime]=FeynDyn(Nbath,finalPoint,deltaKmax,totalT,rho,H,systemCouplingMatrix,w,dw,J,temperature,wholeDensityMatrixOrJustDiagonals,allPointsORjustFinalPoint,cpuORgpu) |
10 | 10 |
|
11 | 11 | %% 1. Fundamental constansts |
12 | 12 | kb=1.3806504*10^(-23); % Joules / Kelvin |
13 | 13 | hbar=1.054571628*10^(-34); % Joules * seconds |
14 | 14 | beta=1/(kb*temperature); |
15 | 15 | %% 2. Setup arrays |
16 | | -M=length(H);M2=M^2;Svector=eig(systemCouplingMatrix).'; |
| 16 | +M=length(H);M2=M^2; |
| 17 | +Svector=eig(systemCouplingMatrix).'; |
| 18 | +Svectors=diag(ones(Nbath,1)); |
17 | 19 | diagonals=diag(reshape(1:M2,M,M)); |
18 | 20 | upperTriangle=find(triu(reshape(1:M2,M,M))); |
19 | 21 | initialPoint=0; |
|
55 | 57 | K=repmat(K,1,length(gridOfTimeIndices)); % for time dependent OQS hamiltonians |
56 | 58 |
|
57 | 59 | [ia,ib] = ndgrid(1:M); |
| 60 | + |
58 | 61 | Sdiff=(kron(Svector,ones(1,M))-kron(ones(1,M),Svector))'; % would not need to transpose the things in the next three lines if Mvector was already a column vector (1:M)' |
59 | 62 |
|
60 | 63 | I_00=expand(exp(-Sdiff.*(kernelEnd(1)*Svector(1,ib)-kernelEndConj(1)*Svector(1,ia)).'),[M2,1]); %eq(14) on pg 4601 of tensor prop. I. repeated values are for n=1, changing values are for n=0 |
61 | 64 | I_mk=zeros(M2^2,deltaKmax+1);I_mkEnd=I_mk;I_mkTD=I_mk;I_mkTDend=I_mk; %I_00 means I0(sk), I_mk means I(sk,s minus k) ? |
62 | 65 | I_mk(:,1)=kron(ones(M2,1),exp(-Sdiff.*(kernel(1)*Svector(1,ib)-kernelConj(1)*Svector(1,ia)).')); %eq(14) on pg 4601 of tensor prop. I. repeated values are for 0 (or just n-1?), changing values are for n %might be better to rearrange stuff so that you can use kron if repmat is slow in fortran, %is it better computationally to do this in one line or separate calcs ? This is delatK=0, ie, eq 14 on pg 4604 of Tensor Prop. I, the diff btwn this and the line above is that this is for sk>=1 while line above is for sk=0 ? eta+kk ? |
63 | 66 | I_mkEnd(:,1)=kron(ones(M2,1),exp(-Sdiff.*(kernelEnd(1)*Svector(1,ib)-kernelEndConj(1)*Svector(1,ia)).')); %repeated values are for 0, changing values are for n %might be better to rearrange stuff so that you can use kron if repmat is slow in fortran, %is it better computationally to do this in one line or separate calcs ? This is delatK=0, eta_NN ? |
| 67 | + |
64 | 68 | for deltaK=1:deltaKmax %should go up to deltaKmax-1 and deltaKmax should be treated separately |
65 | 69 | I_mk(:,1+deltaK)=exp(-kron((kernel(1+deltaK)*Svector(1,ib)-kernelConj(1+deltaK)*Svector(1,ia)).',Sdiff)); %can't we just make kernel a vector and use kernel(deltaK) ? , %saving kernelConj reduces number of times 'conj.m' has to be called, but adds another variable, which is worse ? %is it bad to have such a big argument in kron ? or should i save it and make 2 lines, depending on your units for kernel, you might need to factor the argument by 1/hbar. The reshaping afterwards must be redundant, we must be able to use kron in such a way that the answer ends up in matrix form , |
66 | 70 | I_mkTD(:,1+deltaK)=exp(-kron((kernelTD(1+deltaK)*Svector(1,ib)-kernelTDConj(1+deltaK)*Svector(1,ia)).',Sdiff)); %can't we just make kernel a vector and use kernel(deltaK) ? , %saving kernelConj reduces number of times 'conj.m' has to be called, but adds another variable, which is worse ? %is it bad to have such a big argument in kron ? or should i save it and make 2 lines, depending on your units for kernel, you might need to factor the argument by 1/hbar. The reshaping afterwards must be redundant, we must be able to use kron in such a way that the answer ends up in matrix form , |
67 | 71 | I_mkEnd(:,1+deltaK)=exp(-kron((kernelEnd(1+deltaK)*Svector(1,ib)-kernelEndConj(1+deltaK)*Svector(1,ia)).',Sdiff)); %can't we just make kernel a vector and use kernel(deltaK) ? , %saving kernelConj reduces number of times 'conj.m' has to be called, but adds another variable, which is worse ? %is it bad to have such a big argument in kron ? or should i save it and make 2 lines, depending on your units for kernel, you might need to factor the argument by 1/hbar. The reshaping afterwards must be redundant, we must be able to use kron in such a way that the answer ends up in matrix form , |
68 | 72 | I_mkTDend(:,1+deltaK)=exp(-kron((kernelTDend(1+deltaK)*Svector(1,ib)-kernelTDendConj(1+deltaK)*Svector(1,ia)).',Sdiff)); %can't we just make kernel a vector and use kernel(deltaK) ? , %saving kernelConj reduces number of times 'conj.m' has to be called, but adds another variable, which is worse ? %is it bad to have such a big argument in kron ? or should i save it and make 2 lines, depending on your units for kernel, you might need to factor the argument by 1/hbar. The reshaping afterwards must be redundant, we must be able to use kron in such a way that the answer ends up in matrix form , |
69 | 73 | end %kernel(1)=eta_kk, kernel(2)=eta_{k+1,k} |
| 74 | + |
| 75 | +%% 4.1 Calculation of I tensors with N identical baths |
| 76 | +for jj=2:Nbath |
| 77 | +Svector2=Svectors(Nbath+1-jj,:); %relationship between Svector and Svector2 |
| 78 | +Sdiff2=(kron(Svector2,ones(1,M))-kron(ones(1,M),Svector2))'; |
| 79 | + |
| 80 | +Ij_00=expand(exp(-Sdiff2.*(kernelEnd(1)*Svector2(1,ib)-kernelEndConj(1)*Svector2(1,ia)).'),[M2,1]); |
| 81 | +Ij_mk=zeros(M2^2,deltaKmax+1);Ij_mkEnd=Ij_mk;Ij_mkTD=Ij_mk;Ij_mkTDend=Ij_mk; |
| 82 | +Ij_mk(:,1)=kron(ones(M2,1),exp(-Sdiff2.*(kernel(1)*Svector2(1,ib)-kernelConj(1)*Svector2(1,ia)).')); |
| 83 | +Ij_mkEnd(:,1)=kron(ones(M2,1),exp(-Sdiff2.*(kernelEnd(1)*Svector2(1,ib)-kernelEndConj(1)*Svector2(1,ia)).')); |
| 84 | + |
| 85 | +for deltaK=1:deltaKmax |
| 86 | + Ij_mk(:,1+deltaK)=exp(-kron((kernel(1+deltaK)*Svector2(1,ib)-kernelConj(1+deltaK)*Svector2(1,ia)).',Sdiff2)); |
| 87 | + Ij_mkTD(:,1+deltaK)=exp(-kron((kernelTD(1+deltaK)*Svector2(1,ib)-kernelTDConj(1+deltaK)*Svector2(1,ia)).',Sdiff2)); |
| 88 | + Ij_mkEnd(:,1+deltaK)=exp(-kron((kernelEnd(1+deltaK)*Svector2(1,ib)-kernelEndConj(1+deltaK)*Svector2(1,ia)).',Sdiff2)); |
| 89 | + Ij_mkTDend(:,1+deltaK)=exp(-kron((kernelTDend(1+deltaK)*Svector2(1,ib)-kernelTDendConj(1+deltaK)*Svector2(1,ia)).',Sdiff2)); |
| 90 | +end %kernel(1)=eta_kk, kernel(2)=eta_{k+1,k} |
| 91 | + |
| 92 | +I_00=I_00.*Ij_00; |
| 93 | +I_mk=I_mk.*Ij_mk; |
| 94 | +I_mkTD=I_mkTD.*Ij_mkTD; |
| 95 | +I_mkEnd(:,1)=I_mkEnd(:,1).*Ij_mkEnd(:,1); |
| 96 | +I_mkTDend=I_mkTDend.*Ij_mkTDend; |
| 97 | +end |
| 98 | + |
70 | 99 | %% 5. Propagation of rho for the first deltaKmax timesteps |
71 | 100 | A=K(:,1).*I_mkTD(:,1+1).*I_mk(:,1).*I_00.*expand(rho(:,1),[M2,1]); %I_mkTD=I_k0,k=1 I_mk=I_kk, k=1 , I_00=I_NN,k=0 |
72 | | - |
73 | 101 | Aend=K(:,1).*I_mkTDend(:,1+1).*I_mkEnd(:,1).*I_00.*expand(rho(:,1),[M2,1]); %I_mkTD=I_N0,N=1 I_mkEnd=I_NN, k=1 , I_00=I_NN,N=0 |
74 | | -%K(:,1)=[]; %why is this here ? |
| 102 | +%K(:,1)=[]; % only for time-dependent system Hamiltonians |
75 | 103 |
|
76 | 104 | rho(:,2)=sum(reshape(Aend,M2,M2).'); Aend=[]; |
77 | 105 |
|
78 | | -indices=uint64(npermutek(cast(1:M2,'single'),1+1)); %makes 1:4 single precision. One can make them int(8) but then you'd have to use Jan simons's thing (or you could use Matt Fig's MEX, since you're storing it anyway) |
| 106 | +indices=uint32(npermutek(cast(1:M2,'single'),1+1)); %makes 1:4 single precision. One can make them int(8) but then you'd have to use Jan simons's thing (or you could use Matt Fig's MEX, since you're storing it anyway) |
79 | 107 |
|
80 | | -for J=2:deltaKmax %could go to deltaKmax-1 and incorporate the deltaKmax case into the next forloop, but the former way uses I_mkTD and the latter does not. |
| 108 | +for J=2:deltaKmax %could go to deltaKmax-1 and incorporate the deltaKmax case into the next forloop, but the former way uses I_mkTD and the latter does not. |
81 | 109 |
|
82 | 110 | indices=horzcat(expand(indices,[M2,1]),repmat((1:M2)',size(indices,1),1));% Making 1:M single precision might help |
83 | 111 | A=(expand(A,[M2,1])); |
|
91 | 119 | end |
92 | 120 | A=A.*I_mkTD((indices(:,end-J)-1)*M2+indices(:,end),J+1); %the k=J case for forloop above. uses eta_k0 |
93 | 121 | Aend=Aend.*I_mkTDend((indices(:,end-J)-1)*M2+indices(:,end),J+1); %the k=J case for forloop above. uses eta_N0 |
94 | | - %whos('A','Aend','indices') |
95 | | - |
| 122 | + %whos('A','Aend','indices') |
| 123 | + |
96 | 124 | %weakPaths=find(abs(A)<tol); %in 2001 Sim, it's <= |
97 | 125 | %A(weakPaths)=0;indices(weakPaths,:)=0; % might be better to just remove these rows completely, and then calculate rho using information in the leftover INDICES array. If you remove rows of A, you have to also remove rows of indices to allow A.*I_mk(indices) to work |
98 | 126 | %A=sparse(A);indices=sparse(indices); %now there's a problem with indices=horzcat.... above. Not only does size(indices,1) have to be changed into length(find(indices(:,1)), but expand will also expand the 0's even though you don't want that. |
|
108 | 136 | if strcmpi(cpuORgpu,'GPU');gpuIndex=gpuDevice(1);end |
109 | 137 | KtimesIpermanent=K((indices(:,end-1)-1)*M2+indices(:,end),1); %predefine ? |
110 | 138 | KtimesIpermanentEnd=KtimesIpermanent; |
111 | | -%whos |
| 139 | +%whos |
112 | 140 | for k=0:deltaKmax % it's probably faster if this tensor is just built as A is in the above forloop ... in fact it seems reading indices is slow, so perhaps even the previous iterations should be done like the tensor propagator below: |
113 | 141 | KtimesIpermanent=KtimesIpermanent.*I_mk((indices(:,end-k)-1)*M2+indices(:,end),k+1); %when k=deltaKmax, we're using I_mk(:,1+deltaKmax) which used kernel(1+deltaK) which uses eta 5 |
114 | 142 | KtimesIpermanentEnd=KtimesIpermanentEnd.*I_mkEnd((indices(:,end-k)-1)*M2+indices(:,end),k+1);%when k=deltaKmax, we're using I_mkEnd(:,1+deltaKmax) which used kernelEnd(1+deltaK) which uses eta 3 |
|
127 | 155 | A=reshape(bsxfun(@times,reshape(sum(reshape(A,[],M2),2),1,[]),KtimesIpermanent),[],1); %could be sped up by using some of the calculations in line above |
128 | 156 |
|
129 | 157 | Aend=reshape(Aend,M2,[]); % for some reason, switching the M2 and [] and removing the ,2 in the line below alters the beginning |
130 | | - if or(strcmp(allPointsORjustFinalPoint,'allPoints'),and(strcmp(allPointsORjustFinalPoint,'justFinalPoint'),J==finalPoint)); |
| 158 | + if or(strcmpi(allPointsORjustFinalPoint,'allPoints'),and(strcmpi(allPointsORjustFinalPoint,'justFinalPoint'),J==finalPoint)); |
131 | 159 | switch wholeDensityMatrixOrJustDiagonals |
132 | 160 | case 'justDiagonals' |
133 | 161 | rho(diagonals(1:M-1),J+1)=sum(Aend(diagonals(1:M-1),:),2); % all diagonals but the last |
|
140 | 168 | end |
141 | 169 | end |
142 | 170 | end |
143 | | - |
144 | 171 | elapsedTime=toc; |
145 | 172 | if strcmpi(cpuORgpu,'GPU') |
146 | 173 | rho=gather(rho);reset(gpuIndex); |
|
150 | 177 | function B = expand(A,S) % Compact version of: http://www.mathworks.co.uk/matlabcentral/fileexchange/24536-expand |
151 | 178 |
|
152 | 179 | SA = size(A); % Get the size (and number of dimensions) of input. |
153 | | -if length(SA)~=length(S);error('Length of size vector must equal ndims(A). See help.'); |
| 180 | +if length(SA)~=length(S);error('Length of size vector must equal ndims(A). See help.') |
154 | 181 | elseif any(S~=floor(S));error('The size vector must contain integers only. See help.') |
155 | 182 | end |
156 | 183 | for ii = length(SA):-1:1 |
|
0 commit comments