OpenMCTDHB v2.3

density.F

Go to the documentation of this file.
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 
 All Namespaces Files Functions Variables