OpenMCTDHB v2.3
|
00001 C####################################################################### 00002 C 00003 C MCTDH-module DENSITY 00004 C 00005 C Calculates the reduced density matrices, their inverse matrices and 00006 C their eigenvalues and eigenvectors. 00007 C 00008 C density: passing routine 00009 C denmat: calculates the reduced density matrices 00010 C matinv: inverts a matrix 00011 C 00012 C####################################################################### 00013 00014 C####################################################################### 00015 C DENSITY 00016 C 00017 C passes wavefunction coefficients and relevant part of density matix 00018 C arrays for calculation of reduced density matrices for each mode. 00019 C 00020 C all passed variables are standard arrays, see documentation 00021 C 00022 C V6.3 MB 00023 C####################################################################### 00024 00027 subroutine density(Inv_rho_jk,rho_jk,dicht3,dicht4,workr) 00028 USE moduleinputvariables,ONLY:Morb 00029 implicit none 00030 real*8 dicht3(Morb,2),epsrel,norm,trace,workr(3*MOrb) 00031 REAL*8, PARAMETER :: epsmat = 1.D-8 !KS 00032 complex*16 Inv_rho_jk(Morb,Morb),rho_jk(Morb,Morb), 00033 + dicht4(Morb,Morb) 00034 00035 00036 C----------------------------------------------------------------------- 00037 C compute norm and take a relative regularization parameter 00038 C----------------------------------------------------------------------- 00039 norm=0.0d0 00040 call trhxz(rho_jk,trace,Morb) 00041 norm=norm+trace 00042 epsrel=norm*epsmat 00043 00044 C----------------------------------------------------------------------- 00045 C regularize and invert density matrices for each mode and state 00046 C----------------------------------------------------------------------- 00047 call denmat(Inv_rho_jk,rho_jk, 00048 + dicht3,dicht4,Morb, 00049 + workr,epsrel) 00050 00051 return 00052 end 00053 00054 C####################################################################### 00055 C DENMAT 00056 C 00057 C called by density 00058 C calls to matinv 00059 C 00060 C Calculates the reduced density matrix etc. for a single mode. For 00061 C natural orbital propagation, the eigenvalues and inverse 00062 C density matrices are taken directly from the diagonal of dicht2, 00063 C otherwise dicht2 is diagonalised. 00064 C 00065 C workc: complex workspace needed by subroutine matinv 00066 C for the mode passed by subroutine density: 00067 C dicht1: inverse density matrix 00068 C dicht2: reduced density matrix 00069 C dicht3: eigenvalues and regularised eigenvalues of dicht2 00070 C dicht4: eigenvectors of dicht2 00071 C dim: number of single particle functions 00072 C 00073 C KS 2009 adopted and modified from MCTDH 00074 C####################################################################### 00075 00078 subroutine denmat(dicht1,dicht2,dicht3,dicht4,dim, 00079 + workr,epsrel) 00080 00081 implicit none 00082 00083 integer e,e1,dim 00084 complex*16 dicht1(dim,dim),dicht2(dim,dim), 00085 + dicht4(dim,dim) 00086 real*8 dicht3(dim,2),epsrel,workr(3*dim) 00087 00088 C----------------------------------------------------------------------- 00089 C calculate regularized invers density matrix 00090 C Notice that negative of dicht2 is inverted, so that eigenvalues are 00091 C ordered with largest at dicht3(1) 00092 C----------------------------------------------------------------------- 00093 00094 do e=1,dim 00095 do e1=1,dim 00096 dicht4(e1,e)=-dicht2(e1,e) 00097 enddo 00098 enddo 00099 call matinv(dicht1,dicht3,dicht4,dim,epsrel,workr) 00100 00101 return 00102 end 00103 00104 C####################################################################### 00105 C 00106 C MATINV 00107 C 00108 C called by denmat 00109 C calls to zheev 00110 C 00111 C calls diagonalisation routines for density matrix and from the 00112 C eigenvalues (regularised) and eigenvectors so obtained calculates 00113 C the inverse density matrix 00114 C 00115 C invers: returned as inverted matrix, used as scratch array. 00116 C eval: eigenvalues of passed matrix 00117 C evek: eigenvectors of passed matrix; on input: Matrix to be inverted. 00118 C dim: dimensions of matrices e.g. invers(dim,dim) 00119 C work: complex workspace 00120 C 00121 C KS adopted and modified from MCTDH 00122 C####################################################################### 00123 00126 subroutine matinv(invers,eval,evek,dim,epsrel,work) 00127 00128 implicit none 00129 00130 00131 integer error,e,e1,e2,dim 00132 complex*16 invers(dim,dim),evek(dim,dim) 00133 real*8 eval(dim,2),epsrel,work(3*dim) 00134 00135 C----------------------------------------------------------------------- 00136 C diagonalize matrix 00137 C----------------------------------------------------------------------- 00138 e2 = dim*dim 00139 call zheev('V','U',dim,evek,dim,eval,invers,e2,work,error) 00140 if (error .ne. 0) then 00141 ! routine = 'matinv' 00142 ! message = 'Cannot diagonalize density matrix' 00143 write(*,*)'matinv' 00144 write(*,*)'Cannot diagonalize density matrix' 00145 ! call errormsg 00146 endif 00147 00148 C----------------------------------------------------------------------- 00149 C take negative of eigenvalue (due to re-ordering) and regularize. 00150 C The if-statement is introduced to avoid an underflow. 00151 C----------------------------------------------------------------------- 00152 do e=1,dim 00153 eval(e,1)=-eval(e,1) 00154 if( eval(e,1) .lt. 64.d0*epsrel) then 00155 eval(e,2)=eval(e,1)+epsrel*exp(-eval(e,1)/epsrel) 00156 else 00157 eval(e,2)=eval(e,1) 00158 endif 00159 enddo 00160 00161 C----------------------------------------------------------------------- 00162 C calculate regularized inverse matrix 00163 C----------------------------------------------------------------------- 00164 call zeromxz(invers,dim,dim) 00165 do e=1,dim 00166 do e1=1,dim 00167 do e2=1,dim 00168 invers(e1,e)=invers(e1,e)+ 00169 + dconjg(evek(e,e2))*evek(e1,e2)/eval(e2,2) 00170 enddo 00171 enddo 00172 invers(e,e)=dble(invers(e,e)) 00173 enddo 00174 00175 return 00176 end 00177 00178