|
OpenMCTDHB v2.3
|
|
OpenMCTDHB v2.3
|
00001 C ********************************************************************** 00002 C * * 00003 C * SHORT ITERATIVE LANCZOS (sillib.f) * 00004 C * * 00005 C * Library module containing a short iterative Lanczos integrator. * 00006 C * * 00007 C * Contains: * 00008 C * SILStep: The (real) Lanczos integration routine. * 00009 C * SILErrorMsg: Returns for a given error number a corresponding * 00010 C * error message. * 00011 C * SILDiag: The (real) Lanczos diagonalisation routine. * 00012 C * * 00013 C ********************************************************************** 00014 00015 00016 C ********************************************************************** 00017 C * * 00018 C * SUBROUTINE SILSTEP * 00019 C * * 00020 C * Integrates a system of complex linear first order differential * 00021 C * equations with constant and Hermitian Hamiltonian employing the * 00022 C * short iterative Lanczos method. The routine runs with both * 00023 C * variable step size and order. (First it increases the order to * 00024 C * achieve the desired accuracy. If this failes within the maximum * 00025 C * order, the stepsize is reduced.) SILStep makes only one single * 00026 C * integration step, so it has to be imbedded into a loop that calls * 00027 C * SILStep until the desired time interval is integrated. The ODE is * 00028 C * of the form i dPsi/dt = H|Psi> = Func(Time,Psi,DtPsi) =: DtPsi or * 00029 C * dPsi/dt = -i H|Psi> = Func(Time,Psi,DtPsi) =: DtPsi, depending on * 00030 C * the flag "StdForm". All computations are performed with double * 00031 C * precision. * 00032 C * * 00033 C * Input parameters: * 00034 C * Psi: The (complex) initial-value vector. * 00035 C * DtPsi: Action of Hamiltonian on the initial-value vector, * 00036 C * i.e. H|Psi> or -i H|Psi>, depending on "StdForm". * 00037 C * PsiDim Length of Psi and DtPsi vectors. * 00038 C * IntPeriod: Lenght of time interval to be integrated. * 00039 C * IntOrder: Maximum integration order. * 00040 C * TolError: Maximum error that is tolerated. * 00041 C * Relax: Flag for relaxation calculation. If true, Psi is * 00042 C * relaxated, else propagated. * 00043 C * Restart: Flag for restarting the integrator. If true, the * 00044 C * Krylov space is built up before propagation, else the * 00045 C * old Krylov vectors are used. * 00046 C * StdForm: If true, Func(Time,Psi,DtPsi) = -i H|Psi>, else * 00047 C * Func(Time,Psi,DtPsi) = H|Psi>. * 00048 C * Steps: Number of steps made so far (Steps is passed to * 00049 C * "WriteStep"). * 00050 C * Krylov: Matrix with minimum size PsiDim*(IntOrder-1) with * 00051 C * columns containing the Krylov vectors H^2|psi> to * 00052 C * H^TrueOrder|psi>. The Krylov vectors are needed on * 00053 C * entry only when Restart is false. * 00054 C * * 00055 C * Output parameters: * 00056 C * Psi: Propagated Psi. * 00057 C * DtPsi: Normalised first Krylov vector, i. e. (H-<H>)|Psi>. * 00058 C * Krylov: Matrix with minimum size PsiDim*(IntOrder-1) with * 00059 C * columns containing the Krylov vectors H^2|psi> to * 00060 C * H^TrueOrder|psi>. * 00061 C * Stepsize: Time interval that actually has been integrated (can * 00062 C * be lower than IntPeriod). * 00063 C * TrueOrder: The order that has actually been used (may be less * 00064 C * than IntOrder). * 00065 C * Steps: Same as on entry. * 00066 C * ErrorCode: Error code having the following meaning: * 00067 C * 0: everything was o. k., * 00068 C * 1: illegal integration order, * 00069 C * 2: stepsize underflow, * 00070 C * 3: diagonalization failed, * 00071 C * 4: stepsize too large for restart. * 00072 C * * 00073 C * Other parameters: * 00074 C * CData: Complex*16 data needed in the Func routine. * 00075 C * RData: Real*8 data needed in the Func routine. * 00076 C * IData: Integer data needed in the Func routine. * 00077 C * LData: Logical data needed in the Func routine. * 00078 C * * 00079 C * External routines: * 00080 C * Func: Computes the action of the Hamiltonian on Psi. * 00081 C * Called as * 00082 C * Func(AbsTime,Psi,DtPsi,CData,RData,IData,LData). * 00083 C * WriteStep: Writes the stepsize and error to a file. * 00084 C * Called as WriteStep(Steps,Order,Stepsize,Error). * 00085 C * * 00086 C * Some new Variables: * 00087 C * rlaxnum : This is the argument of the keyword relaxation. * 00088 C * If relaxation has no argument, rlaxnum=-1 . * 00089 C * Relax : logical set true for a relaxation run. * 00090 C * * 00091 C * V6.0 MB * 00092 C * V7.0 GW addition of CData,RData,IData,LData arrays * 00093 C * V8.2 05/01 HDM Addition of 'relaxation' to exited states. * 00094 C * Eigenvector-space now dynamically allocated. * 00095 C ********************************************************************** 00096 00099 Subroutine SILStep(Psi,DtPsi,PsiDim,IntPeriod,IntOrder,TolError, 00100 + Relax,Restart,StdForm,Steps,Krylov,Stepsize, 00101 + TrueOrder,ErrorCode,Time,Func, 00102 + EigenVector,EigenVal,Diagonal,OffDg2,OffDiag, 00103 + CData,RData,IData,LData) 00104 Implicit None 00105 00106 Real*8 RelativeMinStep,Tiny 00107 Parameter (RelativeMinStep = 1.0D-10,Tiny = 1.0D-18) 00108 00109 Logical Relax,Restart,StdForm,LData(*) 00110 Integer PsiDim,IntOrder,TrueOrder,Steps,ErrorCode,IData(*) 00111 Real*8 IntPeriod,TolError,Stepsize,Time,beta1,RData(*) 00112 Complex*16 Psi(PsiDim),DtPsi(PsiDim),Krylov(PsiDim,2:IntOrder), 00113 + CData(*) 00114 External Func 00115 00116 00117 Integer OldOrder,D,P,Q 00118 Real*8 Beta,Error,OldError,MinStepsize,OldStepsize, 00119 + EigenVector(0:IntOrder,0:IntOrder+2), 00120 + Diagonal(0:IntOrder),OffDiag(IntOrder), 00121 + EigenVal(0:IntOrder),OffDg2(IntOrder+1) 00122 Complex*16 Alpha,CBeta,Sum,PreFactor,CInverse 00123 00124 Save OldError,OldStepsize,OldOrder 00125 00126 00127 C --- CHECK INTEGRATION ORDER --- 00128 00129 If ( IntOrder .LT. 2 ) Then 00130 ErrorCode = 1 00131 write(6,'(a,i4,a)') 'IntOrder =',IntOrder,' is too small!' 00132 Return 00133 EndIf 00134 00135 C --- INITIALIZE VARIABLES --- 00136 00137 Stepsize = IntPeriod 00138 MinStepsize = Abs(RelativeMinStep*IntPeriod) 00139 ErrorCode = 0 00140 If (.not.StdForm) Then 00141 PreFactor = (1.0D0,0.0D0) 00142 ElseIf (Relax) Then 00143 PreFactor = (-1.0D0,0.0D0) 00144 Else 00145 PreFactor = (0.0D0,-1.0D0) 00146 EndIf 00147 00148 C --- RETURN RESULT IF DIFFERENTIAL EQUATION IS NOT A SYSTEM --- 00149 00150 If (PsiDim .Eq. 1) Then 00151 TrueOrder = 1 00152 Error = 0.0D0 00153 ! Call WriteStep(Steps,TrueOrder,Stepsize,Error,Time) 00154 Psi(1) = Exp(StepSize*DtPsi(1)/Psi(1))*Psi(1) 00155 Return 00156 EndIf 00157 00158 C --- SKIP CALCULATION OF KRYLOV VECTORS IF DESIRED --- 00159 00160 If (Restart) Then 00161 OldError=0.0D0 00162 OldOrder=0 00163 OldStepsize=0.0D0 00164 Else 00165 If (Abs(Stepsize) .GT. Abs(OldStepsize)+MinStepsize) Then 00166 ErrorCode = 4 00167 Return 00168 EndIf 00169 TrueOrder = OldOrder 00170 Error = Abs(OldError*(Stepsize/OldStepsize)**TrueOrder) 00171 ! Call WriteStep(Steps,TrueOrder,Stepsize,Error,Time) 00172 Goto 200 00173 EndIf 00174 00175 C --- FURTHER INITIALIZE VARIABLES --- 00176 00177 Alpha = (0.0D0,0.0D0) 00178 Beta = 0.0D0 00179 TrueOrder = 1 00180 00181 C --- COMPUTE FIRST DIAGONAL ELEMENT --- 00182 00183 Do D = 1,PsiDim 00184 Alpha = Alpha+DConjg(Psi(D))*DtPsi(D) 00185 EndDo 00186 Diagonal(0) = Alpha/PreFactor 00187 00188 C --- DETERMINE CORRESPONDING BASIS VECTOR AND OFF-DIAGONAL ELEMENT --- 00189 00190 Do D = 1,PsiDim 00191 DtPsi(D) = DtPsi(D)-Alpha*Psi(D) 00192 Beta = Beta+Dble(DConjg(DtPsi(D))*DtPsi(D)) 00193 EndDo 00194 00195 Beta = Sqrt(Beta) 00196 OffDiag(1) = Beta 00197 beta1 = Beta 00198 CBeta = PreFactor*Beta 00199 00200 C --- NORMALIZE BASIS VECTOR --- 00201 00202 If (Beta .LT. 1.d-20) Then 00203 Beta = 0.0D0 00204 CInverse = (0.0D0,0.0D0) 00205 Else 00206 CInverse = 1.0D0/CBeta 00207 EndIf 00208 Do D = 1,PsiDim 00209 DtPsi(D) = CInverse*DtPsi(D) 00210 EndDo 00211 00212 C --- COMPUTE ERROR ESTIMATE --- 00213 00214 Error = Abs(Beta*Stepsize) 00215 00216 C --- WRITE STEPSIZE AND ERROR IF DESIRED --- 00217 00218 ! Call WriteStep(Steps,TrueOrder,Stepsize,Error,Time) 00219 00220 C --- BUILD UP THE KRYLOV SPACE --- 00221 00222 100 Continue 00223 00224 C --- CHECK IF STEPSIZE IS TOO SMALL --- 00225 00226 If (Stepsize .LT. MinStepSize) Then 00227 ErrorCode = 2 00228 Return 00229 EndIf 00230 00231 C --- RE-INITIALIZE VARIABLES --- 00232 00233 Alpha = (0.0D0,0.0D0) 00234 Beta = 0.0D0 00235 TrueOrder = TrueOrder+1 00236 00237 00238 C --- EVALUATE FUNCTION WITH LAST BASIS VECTOR --- 00239 00240 If (TrueOrder .Eq. 2) Then 00241 Call Func(Time,DtPsi,Krylov(1,2),CData,RData,IData,LData) 00242 Else 00243 Call Func(Time,Krylov(1,TrueOrder-1),Krylov(1,TrueOrder), 00244 + CData,RData,IData,LData) 00245 EndIf 00246 00247 C --- COMPUTE DIAGONAL ELEMENT --- 00248 00249 C Note that the Krylov vectors number zero and one aren't stored in 00250 C "Krylov" but in "Psi" and "DtPsi", respectively. 00251 00252 If (TrueOrder .Eq. 2) Then 00253 Do D = 1,PsiDim 00254 Alpha = Alpha+DConjg(DtPsi(D))*Krylov(D,2) 00255 EndDo 00256 Else 00257 Do D = 1,PsiDim 00258 Alpha = Alpha+DConjg(Krylov(D,TrueOrder-1)) 00259 + *Krylov(D,TrueOrder) 00260 EndDo 00261 EndIf 00262 Diagonal(TrueOrder-1) = Alpha/PreFactor 00263 00264 C --- COMPUTE OFF-DIAGONAL ELEMENT AND BASIS VECTOR --- 00265 00266 If (TrueOrder .Eq. 2) Then 00267 Do D = 1,PsiDim 00268 Krylov(D,2) = Krylov(D,2)-Alpha*DtPsi(D)-CBeta*Psi(D) 00269 Beta = Beta+Dble(DConjg(Krylov(D,2))*Krylov(D,2)) 00270 EndDo 00271 ElseIf (TrueOrder .Eq. 3) Then 00272 Do D = 1,PsiDim 00273 Krylov(D,3) = Krylov(D,3)-Alpha*Krylov(D,2)-CBeta*DtPsi(D) 00274 Beta = Beta+Dble(DConjg(Krylov(D,3))*Krylov(D,3)) 00275 EndDo 00276 Else 00277 Do D = 1,PsiDim 00278 Krylov(D,TrueOrder) = Krylov(D,TrueOrder) 00279 + -Alpha*Krylov(D,TrueOrder-1)-CBeta*Krylov(D,TrueOrder-2) 00280 Beta = Beta+DConjg(Krylov(D,TrueOrder))*Krylov(D,TrueOrder) 00281 EndDo 00282 EndIf 00283 00284 Beta = Sqrt(Beta) 00285 OffDiag(TrueOrder) = Beta 00286 CBeta = PreFactor*Beta 00287 00288 C --- NORMALIZE BASIS VECTOR --- 00289 00290 If (Beta .LT. 1.d-20) Then 00291 OffDiag(TrueOrder) = 0.0D0 00292 Beta = 0.0D0 00293 CInverse = (0.0D0,0.0D0) 00294 Else 00295 CInverse = 1.0D0/CBeta 00296 EndIf 00297 00298 Do D = 1,PsiDim 00299 Krylov(D,TrueOrder) = CInverse*Krylov(D,TrueOrder) 00300 EndDo 00301 00302 C --- COMPUTE ERROR ESTIMATE --- 00303 00304 Error = Abs(Error*Beta*Stepsize/Dble(TrueOrder)) 00305 00306 C --- WRITE STEPSIZE AND ERROR IF DESIRED --- 00307 00308 ! Call WriteStep(Steps,TrueOrder,Stepsize,Error,Time) 00309 00310 C --- CONTINUE IF NOT CONVERGED --- 00311 00312 If((TrueOrder.LT.IntOrder) .And. (Error.GT.TolError) 00313 + .And. (Beta.gt.1.d-20) ) Goto 100 00314 00315 C --- ESTIMATE NEXT MATRIX ELEMENT --- 00316 00317 Diagonal(TrueOrder) = Diagonal(TrueOrder-1) 00318 c Diagonal(TrueOrder) = Diagonal(TrueOrder-1) + beta 00319 00320 C --- DECREASE STEPSIZE IF LANCZOS DIDN'T CONVERGE --- 00321 00322 If (Error .GT. TolError) Then 00323 Stepsize = Stepsize*(TolError/Error)**(1.0D0/Dble(TrueOrder)) 00324 ! Call WriteStep(Steps,TrueOrder,Stepsize,TolError,Time) 00325 EndIf 00326 00327 00328 C --- SAVE SOME VALUES IN CASE THAT NEXT CALL IS NON-RESTART CALL --- 00329 00330 OldOrder = TrueOrder 00331 OldStepsize = Stepsize 00332 OldError = Error 00333 00334 C --- Save Diagonal and off diagonal values since the diagonalisation --- 00335 C --- routine will overwrite those. --- 00336 00337 call cpvxd(Diagonal,EigenVal,TrueOrder+1) 00338 call cpvxd(OffDiag,OffDg2,TrueOrder) 00339 00340 C --- DIAGONALIZE LANCZOS MATRIX --- CALL LAPACK ROUTINE --- 00341 C --- EigenVector(*,IntOrder+1) serves as work array 00342 if(beta.lt.1.d-18) TrueOrder=TrueOrder-1 00343 call DSTEQR('I',TrueOrder+1,EigenVal,OffDg2,EigenVector, 00344 + IntOrder+1,EigenVector(1,IntOrder+1),ErrorCode) 00345 If (ErrorCode .NE. 0) Then 00346 write(*,*)"DSTEQR ERRORCODE",ErrorCode 00347 ErrorCode = 3 00348 Return 00349 EndIf 00350 00351 00352 200 Continue ! Jump to here, if not RESTART. 00353 00354 C --- PROPAGATE WAVEFUNCTION --- 00355 00356 C Note again that the Krylov vectors number zero and one aren't stored 00357 C in "Krylov" but in "Psi" and "DtPsi", respectively. 00358 00359 Do P = 0,TrueOrder 00360 Sum = (0.0D0,0.0D0) 00361 Do Q = 0,TrueOrder 00362 Sum = Sum+Exp(PreFactor*EigenVal(Q)*Stepsize)* 00363 + EigenVector(0,Q)*EigenVector(P,Q) 00364 EndDo 00365 If (P .Eq. 0) Then 00366 Do D = 1,PsiDim 00367 Psi(D) = Sum*Psi(D) 00368 EndDo 00369 ElseIf (P .Eq. 1) Then 00370 Do D = 1,PsiDim 00371 Psi(D) = Psi(D)+Sum*DtPsi(D) 00372 EndDo 00373 Else 00374 Do D = 1,PsiDim 00375 Psi(D) = Psi(D)+Sum*Krylov(D,P) 00376 EndDo 00377 EndIf 00378 EndDo 00379 00380 Return 00381 End 00382 00383 00384 C ********************************************************************** 00385 C * * 00386 C * SUBROUTINE SILERRORMSG * 00387 C * * 00388 C * Generates for a given error number returned by "SILStep/SILStep2/ * 00389 C * SILStep3" a corresponding error message. * 00390 C * * 00391 C * Input parameters: * 00392 C * Error: Error code returned by SILStep/SILStep2. * 00393 C * * 00394 C * Output parameters: * 00395 C * Msg: Error message. * 00396 C * * 00397 C * V7.0 MB * 00398 C * * 00399 C ********************************************************************** 00400 00403 Subroutine SILErrorMsg (Error,Msg) 00404 00405 Implicit None 00406 00407 Integer Error 00408 Character*(*) Msg 00409 00410 C --- GENERATE ERROR MESSAGE --- 00411 00412 If (Error .Eq. 1) Then 00413 Msg = 'Illegal integration order' 00414 ElseIf (Error .Eq. 2) Then 00415 Msg = 'Stepsize underflow' 00416 ElseIf (Error .Eq. 3) Then 00417 Msg = 'Diagonalisation failed' 00418 ElseIf (Error .Eq. 4) Then 00419 Msg = 'Illegal initial stepsize' 00420 Else 00421 Msg = 'Unknown error occurred' 00422 EndIf 00423 00424 Return 00425 End 00426 00427 00428 C ********************************************************************** 00429 C * * 00430 C * SUBROUTINE SILDIAG * 00431 C * * 00432 C * Diagonalises a (complex) Hermitian matrix by Lanczos-iterations. * 00433 C * This routine is used for improved relaxation and for meigenf. * 00434 C * * 00435 C * Input parameters: * 00436 C * Psi: The (complex) initial-value vector. * 00437 C * DtPsi: Action of Hamiltonian on the initial-value vector, * 00438 C * i.e. H|Psi> or -i H|Psi>, depending on "StdForm". * 00439 C * PsiDim Length of Psi and DtPsi vectors. * 00440 C * IntPeriod: Lenght of time interval to be integrated. * 00441 C * IntOrder: Maximum integration order. * 00442 C * TolError: Maximum error that is tolerated. * 00443 C * StdForm: If true, Func(Time,Psi,DtPsi) = -i H|Psi>, else * 00444 C * Func(Time,Psi,DtPsi) = H|Psi>. * 00445 C * Steps: Number of steps made so far (Steps is passed to * 00446 C * "WriteStep"). * 00447 C * Krylov: Matrix with minimum size PsiDim*(IntOrder-1) with * 00448 C * columns containing the Krylov vectors H^2|psi> to * 00449 C * H^TrueOrder|psi>. * 00450 C * TrueOrder: On input: the minimal order to be used. * 00451 C * ortho: if .true., then a full re-orthogonalisation of the * 00452 C * Krylov vectors is performed. * 00453 C * * 00454 C * Output parameters: * 00455 C * Psi: Computed eigenvector. The number of the eigenvector * 00456 C * to be chosen is specified by rlaxnum. * 00457 C * DtPsi: Normalised first Krylov vector, i. e. (H-<H>)|Psi>. * 00458 C * Krylov: Matrix with minimum size PsiDim*(IntOrder-1) with * 00459 C * columns containing the Krylov vectors H^2|psi> to * 00460 C * H^TrueOrder|psi>. * 00461 C * TrueOrder: The order that has actually been used (may be less * 00462 C * than IntOrder). * 00463 C * Steps: Same as on entry. * 00464 C * ErrorCode: Error code having the following meaning: * 00465 C * 0: everything was o. k., * 00466 C * 1: illegal integration order, * 00467 C * 3: diagonalization failed, * 00468 C * * 00469 C * Other parameters: * 00470 C * CData: Complex*16 data needed in the Func routine. * 00471 C * RData: Real*8 data needed in the Func routine. * 00472 C * IData: Integer data needed in the Func routine. * 00473 C * LData: Logical data needed in the Func routine. * 00474 C * * 00475 C * External routines: * 00476 C * Func: Computes the action of the Hamiltonian on Psi. * 00477 C * Called as * 00478 C * Func(AbsTime,Psi,DtPsi,CData,RData,IData,LData). * 00479 C * WriteStep: Writes the stepsize and error to a file. * 00480 C * Called as WriteStep(Steps,Order,Stepsize,Error). * 00481 C * * 00482 C * Some new Variables: * 00483 C * rlaxnum : This is the argument of the keyword relaxation. * 00484 C * If relaxation has no argument, rlaxnum=-1 . * 00485 C * * 00486 C * V8.3 10/03 HDM * 00487 C ********************************************************************** 00488 00491 Subroutine SILdiag(Psi,DtPsi,PsiDim,IntOrder,TolError,StdForm, 00492 + Steps,Krylov,TrueOrder,ErrorCode,rlaxnum,Time, 00493 + ortho,Func,EigenVector,EigenVal,Diagonal, 00494 + OffDg2,OffDiag,CData,RData,IData,LData) 00495 Implicit None 00496 00497 00498 Logical StdForm,FirstCall,full,ortho,LData(*) 00499 Integer PsiDim,IntOrder,TrueOrder,Steps,ErrorCode,rlaxnum, 00500 + IData(*) 00501 Real*8 TolError,Stepsize,Time,beta1,RData(*) 00502 Complex*16 Psi(PsiDim),DtPsi(PsiDim),Krylov(PsiDim,2:IntOrder), 00503 + CData(*) 00504 External Func 00505 00506 00507 Integer minorder,D,P,Q,i,nd,relaxnum 00508 Real*8 Beta,Error,mxvc,err1,Tinit, 00509 + EigenVector(0:IntOrder,0:IntOrder+2), 00510 + Diagonal(0:IntOrder),OffDiag(IntOrder), 00511 + EigenVal(0:IntOrder),OffDg2(IntOrder+1) 00512 Complex*16 Alpha,CBeta,PreFactor,CInverse,ovlp 00513 00514 data FirstCall/.true./ 00515 Save FirstCall, full 00516 Save Tinit,relaxnum,minorder 00517 00518 C --- CHECK INTEGRATION ORDER --- 00519 00520 If ( IntOrder .LT. 2 ) Then 00521 ErrorCode = 1 00522 write(6,'(a,i4,a)') 'IntOrder =',IntOrder,' is too small!' 00523 Return 00524 EndIf 00525 00526 C --- INITIALIZE VARIABLES --- 00527 00528 nd = 2 00529 Stepsize = 0.d0 00530 ErrorCode = 0 00531 full = .false. 00532 if(FirstCall) Then 00533 Tinit = Time 00534 00535 if(rlaxnum .ge. 9999) then 00536 full = .true. 00537 relaxnum = rlaxnum - 10000 00538 nd = IntOrder+1 00539 else 00540 relaxnum = rlaxnum 00541 full = .false. 00542 endif 00543 If (relaxnum .lt. 0) Then 00544 Write(6,*) ' relaxnum =',relaxnum, ' is .lt. 0. STOP' 00545 stop 00546 EndIf 00547 If (.not.StdForm) Then ! This is for meigenf 00548 minorder = min(TrueOrder,IntOrder) 00549 Else 00550 minorder = 0 00551 EndIf 00552 endif 00553 If (.not.StdForm) Then 00554 PreFactor = (1.0D0,0.0D0) 00555 Else 00556 PreFactor = (-1.0D0,0.0D0) 00557 EndIf 00558 00559 C --- FURTHER INITIALIZE VARIABLES --- 00560 00561 Alpha = (0.0D0,0.0D0) 00562 Beta = 0.0D0 00563 TrueOrder = 1 00564 00565 C --- COMPUTE FIRST DIAGONAL ELEMENT --- 00566 00567 Do D = 1,PsiDim 00568 Alpha = Alpha+DConjg(Psi(D))*DtPsi(D) 00569 EndDo 00570 Diagonal(0) = Alpha/PreFactor 00571 00572 C --- DETERMINE CORRESPONDING BASIS VECTOR AND OFF-DIAGONAL ELEMENT --- 00573 00574 Do D = 1,PsiDim 00575 DtPsi(D) = DtPsi(D)-Alpha*Psi(D) 00576 Beta = Beta+Dble(DConjg(DtPsi(D))*DtPsi(D)) 00577 EndDo 00578 00579 C --- Re-orthogonalize DtPsi 00580 if(ortho) then 00581 ovlp = (0.d0,0.d0) 00582 Do D = 1,PsiDim 00583 ovlp = ovlp + DConjg(Psi(D))*DtPsi(D) 00584 enddo 00585 Beta = 0.0D0 00586 Do D = 1,PsiDim 00587 DtPsi(D)=DtPsi(D)-ovlp*Psi(D) 00588 Beta = Beta+Dble(DConjg(DtPsi(D))*DtPsi(D)) 00589 enddo 00590 endif 00591 00592 Beta = Sqrt(Beta) 00593 OffDiag(1) = Beta 00594 beta1 = Beta 00595 CBeta = PreFactor*Beta 00596 00597 C --- NORMALIZE BASIS VECTOR --- 00598 00599 If (Beta .LT. 1.d-20) Then 00600 Beta = 0.0D0 00601 CInverse = (0.0D0,0.0D0) 00602 Else 00603 CInverse = 1.0D0/CBeta 00604 EndIf 00605 Do D = 1,PsiDim 00606 DtPsi(D) = CInverse*DtPsi(D) 00607 EndDo 00608 00609 C --- COMPUTE ERROR ESTIMATE --- 00610 00611 Error = Beta 00612 00613 C --- WRITE STEPSIZE AND ERROR IF DESIRED --- 00614 00615 ! Call WriteStep(Steps,TrueOrder,Stepsize,Error,Time) 00616 00617 C --- BUILD UP THE KRYLOV SPACE --- 00618 00619 100 Continue 00620 00621 C --- RE-INITIALIZE VARIABLES --- 00622 00623 Alpha = (0.0D0,0.0D0) 00624 Beta = 0.0D0 00625 TrueOrder = TrueOrder+1 00626 00627 00628 C --- EVALUATE FUNCTION WITH LAST BASIS VECTOR --- 00629 00630 If (TrueOrder .Eq. 2) Then 00631 Call Func(Time,DtPsi,Krylov(1,2),CData,RData,IData,LData) 00632 Else 00633 Call Func(Time,Krylov(1,TrueOrder-1),Krylov(1,TrueOrder), 00634 + CData,RData,IData,LData) 00635 EndIf 00636 00637 C --- COMPUTE DIAGONAL ELEMENT --- 00638 00639 C Note that the Krylov vectors number zero and one aren't stored in 00640 C "Krylov" but in "Psi" and "DtPsi", respectively. 00641 00642 If (TrueOrder .Eq. 2) Then 00643 Do D = 1,PsiDim 00644 Alpha = Alpha+DConjg(DtPsi(D))*Krylov(D,2) 00645 EndDo 00646 Else 00647 Do D = 1,PsiDim 00648 Alpha = Alpha+DConjg(Krylov(D,TrueOrder-1)) 00649 + *Krylov(D,TrueOrder) 00650 EndDo 00651 EndIf 00652 Diagonal(TrueOrder-1) = Alpha/PreFactor 00653 00654 C --- COMPUTE OFF-DIAGONAL ELEMENT AND BASIS VECTOR --- 00655 00656 If (TrueOrder .Eq. 2) Then 00657 Do D = 1,PsiDim 00658 Krylov(D,2) = Krylov(D,2)-Alpha*DtPsi(D)-CBeta*Psi(D) 00659 Beta = Beta+Dble(DConjg(Krylov(D,2))*Krylov(D,2)) 00660 EndDo 00661 ElseIf (TrueOrder .Eq. 3) Then 00662 Do D = 1,PsiDim 00663 Krylov(D,3) = Krylov(D,3)-Alpha*Krylov(D,2)-CBeta*DtPsi(D) 00664 Beta = Beta+Dble(DConjg(Krylov(D,3))*Krylov(D,3)) 00665 EndDo 00666 Else 00667 Do D = 1,PsiDim 00668 Krylov(D,TrueOrder) = Krylov(D,TrueOrder) 00669 + -Alpha*Krylov(D,TrueOrder-1)-CBeta*Krylov(D,TrueOrder-2) 00670 Beta = Beta+DConjg(Krylov(D,TrueOrder))*Krylov(D,TrueOrder) 00671 EndDo 00672 EndIf 00673 00674 00675 C --- Re-Orthogonalize against Psi and DtPsi --- 00676 ovlp =(0.d0,0.d0) 00677 Do D = 1,PsiDim 00678 ovlp = ovlp + DConjg(Psi(D))*Krylov(D,TrueOrder) 00679 enddo 00680 Do D = 1,PsiDim 00681 Krylov(D,TrueOrder)=Krylov(D,TrueOrder)-ovlp*Psi(D) 00682 enddo 00683 ovlp =(0.d0,0.d0) 00684 Do D = 1,PsiDim 00685 ovlp = ovlp + DConjg(DtPsi(D))*Krylov(D,TrueOrder) 00686 enddo 00687 Do D = 1,PsiDim 00688 Krylov(D,TrueOrder)=Krylov(D,TrueOrder)-ovlp*DtPsi(D) 00689 enddo 00690 00691 00692 C --- Re-Orthogonalize against Krylov --- 00693 if((ortho.or.full) .and. TrueOrder .ge. 3) then 00694 Do i = 2, TrueOrder-1 00695 ovlp =(0.d0,0.d0) 00696 Do D = 1,PsiDim 00697 ovlp = ovlp + DConjg(Krylov(D,i))*Krylov(D,TrueOrder) 00698 enddo 00699 Do D = 1,PsiDim 00700 Krylov(D,TrueOrder) = Krylov(D,TrueOrder) - 00701 + ovlp*Krylov(D,i) 00702 enddo 00703 enddo 00704 00705 C --- Determine Beta 00706 Beta = 0.d0 00707 Do D = 1,PsiDim 00708 Beta = Beta+DConjg(Krylov(D,TrueOrder))* 00709 + Krylov(D,TrueOrder) 00710 enddo 00711 endif 00712 00713 Beta = Sqrt(Beta) 00714 OffDiag(TrueOrder) = Beta 00715 CBeta = PreFactor*Beta 00716 00717 C --- NORMALIZE BASIS VECTOR --- 00718 00719 If (Beta .LT. 1.d-20) Then 00720 OffDiag(TrueOrder) = 0.0D0 00721 Beta = 0.0D0 00722 CInverse = (0.0D0,0.0D0) 00723 Else 00724 CInverse = 1.0D0/CBeta 00725 EndIf 00726 00727 Do D = 1,PsiDim 00728 Krylov(D,TrueOrder) = CInverse*Krylov(D,TrueOrder) 00729 EndDo 00730 00731 C --- COMPUTE ERROR ESTIMATE --- 00732 00733 Error = Error*Beta/(0.5d0*Abs(alpha-Diagonal(0))+Beta) 00734 00735 C --- WRITE STEPSIZE AND ERROR IF DESIRED --- 00736 00737 ! Call WriteStep(Steps,TrueOrder,Stepsize,Error,Time) 00738 00739 C --- CONTINUE IF NOT CONVERGED --- 00740 00741 If((TrueOrder.LT.IntOrder) .And. (full.or.(Error.GT.TolError)) 00742 + .and. Beta.gt.1.d-20 .or. TrueOrder.LT.MinOrder ) Goto 100 00743 00744 C --- ESTIMATE NEXT MATRIX ELEMENT --- 00745 00746 c Diagonal(TrueOrder) = Diagonal(TrueOrder-1) 00747 Diagonal(TrueOrder) = Diagonal(TrueOrder-1) + beta 00748 00749 00750 C --- SET nd TO AVOID EXCESSIVE TRIAL DIAGONALISATIONS --- 00751 00752 if(TrueOrder.LT.IntOrder .And. 00753 + Beta.gt.1.d-20 .And. mod(TrueOrder,nd) .ne. 0 ) goto 100 00754 if(TrueOrder.ge. 20) nd = 5 ! This is to avoid excessive 00755 if(TrueOrder.ge.100) nd = 10 ! trial diagonalisations. 00756 if(TrueOrder.ge.200) nd = 20 00757 if(TrueOrder.ge.300) nd = 40 00758 if(TrueOrder.ge.500) nd = 80 00759 00760 C --- Save Diagonal and off diagonal values since the diagonalisation --- 00761 C --- routine will overwrite those. --- 00762 00763 call cpvxd(Diagonal,EigenVal,TrueOrder+1) 00764 call cpvxd(OffDiag,OffDg2,TrueOrder) 00765 00766 C --- DIAGONALIZE LANCZOS MATRIX --- CALL LAPACK ROUTINE --- 00767 C --- EigenVector(*,IntOrder+1) serves as work array 00768 if(beta.lt.1.d-18) TrueOrder=TrueOrder-1 00769 call DSTEQR('I',TrueOrder+1,EigenVal,OffDg2,EigenVector, 00770 + IntOrder+1,EigenVector(1,IntOrder+1),ErrorCode) 00771 00772 If (ErrorCode .NE. 0) Then 00773 write(*,*)"DSTEQR ERRORCODE",ErrorCode 00774 write(*,*)"TRUEORDER",TRUEORDER 00775 write(*,*)"EigenVal",EigenVal 00776 ErrorCode = 3 00777 Return 00778 EndIf 00779 00780 00781 00782 C --- RELAX WAVEFUNCTION TO EIGENSTATE NO. Q --- 00783 C --- SEARCH FOR THE EIGENVECTOR THAT IS CLOSEST TO PSI --- 00784 mxvc = 0.d0 00785 Do P = 0,TrueOrder 00786 OffDg2(p+1) = EigenVector(0,P)**2 00787 If(OffDg2(p+1) .gt. mxvc ) Then 00788 mxvc = OffDg2(p+1) 00789 Q = P 00790 EndIf 00791 EndDo 00792 00793 if(time.eq.tinit .and. relaxnum.le.900) 00794 + Q = min(relaxnum,TrueOrder) 00795 err1 = 1000.d0*Beta*EigenVector(TrueOrder,Q)**2 00796 if(.not.StdForm) ! for meigenf 00797 + err1 = err1 + 100.d0*Beta*EigenVector(TrueOrder,minorder/4)**2 00798 + + 100.d0*Beta*EigenVector(TrueOrder,minorder/2)**2 00799 error = max(error,err1) 00800 00801 if( TrueOrder.LT.IntOrder .And. Beta.gt.1.d-20 00802 + .And. err1.gt.TolError ) goto 100 00803 00804 C --- Build new A-vector. 00805 Do D = 1,PsiDim 00806 Psi(D) = EigenVector(0,Q)*Psi(D)+EigenVector(1,Q)*DtPsi(D) 00807 EndDo 00808 Do P = 2,TrueOrder 00809 Do D = 1,PsiDim 00810 Psi(D) = Psi(D)+EigenVector(P,Q)*Krylov(D,P) 00811 EndDo 00812 EndDo 00813 00814 C --- Write to rlx_info file. EigenVector(*,IntOrder+1) is a work array, 00815 C containing the error-estimates. 00816 if(FirstCall) then 00817 do P = 0,TrueOrder 00818 EigenVector(p,IntOrder+1) = 00819 + 2.d0*27.2114d0*Beta*EigenVector(TrueOrder,p)**2 00820 C........The extra factor of 2 is due to the modification, line 358. 00821 if(EigenVector(p,IntOrder+1).lt.1.d-30) 00822 + EigenVector(p,IntOrder+1) = 0.d0 00823 enddo 00824 else 00825 EigenVector(Q,IntOrder+1) = 00826 + 2.d0*27.2114d0*Beta*EigenVector(TrueOrder,Q)**2 00827 endif 00828 !KS call wrtrlx(Q,TrueOrder,mxvc,beta1,time,EigenVal, 00829 !KS + EigenVector(0,IntOrder+1),OffDg2,psi,FirstCall,StdForm) 00830 if(StdForm) then 00831 FirstCall = .false. ! Normal case, improved relaxation 00832 else 00833 steps = Q ! This is for meigenf 00834 do p =0, IntOrder 00835 EigenVector(p,IntOrder+2)=EigenVal(p) 00836 enddo 00837 endif 00838 00839 Return 00840 End 00841 00842 00845 subroutine WriteStep(Step,Order,Stepsize,Error,Time) 00846 !---{{{ 00847 Implicit None 00848 00849 Real*8 Stepsize,Error,Time 00850 Integer Step,Order 00851 00852 ! --- WRITE STEPSIZE AND ERROR --- 00853 00854 !open(unit=1981,file='steps.dat',form='FORMATTED',& 00855 ! status='UNKNOWN',position='APPEND') 00856 ! 00857 ! Write(1981,'(I7,I13,2E19.6E3,2E19.6E3,F18.3)') & 00858 ! Step, Order, Stepsize, Error,Time 00859 ! 00860 !close(1981) 00861 00862 Return 00863 !---}}} 00864 end subroutine WriteStep 00865 00866 00867
1.7.3
