VASP development version

Here we document our internal development version of vasp. In order for something to be added to this version, there should be a reasonable amount of generality, and there should be plenty of comments and or documentation so others can figure out what you have done. Any modification of the code should start and end with a BEGIN-MODVASP:# and END-MODVASP:# tag, where # is an integer assigned to the give change that you have made. There is a git repository for maintaining this source code, and this is strictly limited to group members. You can clone the repsitory as follows:

$ git clone ssh://grandcentral/home/cam1/git_repos/modvasp.git

The above assumes that grandcentral is set up in your .ssh/config file. The makefile is the pure intel version. The vasp.5.lib directory is unchanged. Eventually we will have a binary that is on all machines in /usr/local/bin/ but I have not had time to do this yet.

Reading in density matrix

When you are running DFT+U, you will often want to initialize in some specific density matrix for the correlated subspace. In order to do this, we have modified LDApU.F in two ways (modifications # 1 and 2). Modification 2 keeps track of how many density matrix construction iterations have taken place. We do this by simply reading/writing a file as follows:

! these are the extra variables that we had to declare...
logical iffile,isreal,proceed
real(q)   tempdm(2*LANG_LDAPLUSU(ITYP)+1,   2*LANG_LDAPLUSU(ITYP)+1)
real(q) tempdm_c(2*LANG_LDAPLUSU(ITYP)+1,2*(2*LANG_LDAPLUSU(ITYP)+1))
integer it,lcor
character (len=90) :: fn,dmfile

! BEGIN-MODVASP:2 - Keep track of how many times this subroutine has been called...
! This is not the total number of electronic iterations as this routine is not called
! every electronic step. We only increment the iteration once for each density update.
!
! If the file does not exist then create it and set iteration to 1...
inquire(file='iteration.out',exist=iffile)
if (iffile.eq. .false.)then
    open(65,file="iteration.out")
    write(65,*) '2'
    close(65)
    it=1
else  ! read the file...
  if (IATOM.eq.1)then ! only increment if IATOM=1...
    open(65,file="iteration.out")
    read(65,*) it
    close(65)
    open(65,file="iteration.out")
    write(65,*) it+1
    close(65)
  else ! just read the iteration number for other atoms...
    open(65,file="iteration.out")
    read(65,*) it
    it=it-1 ! we have to decrement this because IATOM=1 already incremented it...
    close(65)
  endif
endif
! END-MODVASP:2

Then next step is modification #1 where we actually read in the density matrix. You can create an initial iteration.out file, and if you make the integer negative the density matrix will be read in until it becomes positive. There are several options to consider. First, if spin-orbit coupling is used then the density matrix is complex. The other factor to consider is that there may be more than one correlated atom per unit cell. Therefore, there can be four different kinds of input files:

  1. density_matrix.inp - This will be for the real density matrix and it is assumed that all correlated atoms will have the same density matrix. You will need to enter two (ie. one for each spin) 5x5 blocks of real numnbers for d-electrons or two blocks of 7x7 real numbers for f-electrons. To be painfully clear, for the case of f-electrons we would have:

 0.9118  0.0225 -0.1257  0.2177 -0.0857  0.0464  0.1445
 0.0225  0.2372  0.0187  0.0148 -0.0170 -0.0074  0.0089
-0.1257  0.0187  0.8592  0.3208 -0.0114  0.0762  0.0373
 0.2177  0.0148  0.3208  0.3914 -0.0541 -0.1417  0.0573
-0.0857 -0.0170 -0.0114 -0.0541  0.7715 -0.0036  0.3986
 0.0464 -0.0074  0.0762 -0.1417 -0.0036  0.9954  0.0134
 0.1445  0.0089  0.0373  0.0573  0.3986  0.0134  0.4062

 0.0345 -0.0004  0.0074 -0.0010  0.0005 -0.0007 -0.0008
-0.0004  0.1287 -0.0003  0.0003 -0.0007  0.0000  0.0009
 0.0074 -0.0003  0.0317 -0.0016  0.0001  0.0000 -0.0003
-0.0010  0.0003 -0.0016  0.0416 -0.0001  0.0010  0.0001
 0.0005 -0.0007  0.0001 -0.0001  0.0310  0.0001 -0.0079
-0.0007  0.0000  0.0000  0.0010  0.0001  0.0254 -0.0001
-0.0008  0.0009 -0.0003  0.0001 -0.0079 -0.0001  0.0360

  1. density_matrix.inp.XX - This is the same as above except XX will be an integer specifying a specific atom. This is useful if you want to set up an orbital density wave. This file takes precedence over density_matrix.inp.

  2. density_matrix_comp.inp - Same as the first option except each row has to have another block with all of the imaginary entries. This follows the same format that vasp uses to output the complex density matrix.

  3. density_matrix_comp.inp.XX - Same as option 2 except you can specify a given atom.

! BEGIN-MODVASP:1 - read in the density matrix if desired...

! define the correlated subspace: lcor=3 for f, 2 for d, etc.
lcor=LANG_LDAPLUSU(ITYP)

! NCDIJ is the number of spins (ie. 2 for normal sp calc)...
if (it <0)then
  proceed=.FALSE.
  print*,'***Read in density matrix***'

  ! look for atom specific real density matrix file...
  write (fn, '( "density_matrix.inp.", i2.2 )') IATOM
  inquire(file=fn,exist=iffile)
  if (iffile.eq. .true.) then
    isreal=.TRUE.;dmfile=fn;proceed=.TRUE.
  endif

  ! look for atom specific complex density matrix file...
  write (fn, '( "density_matrix_comp.inp.", i2.2 )') IATOM
  inquire(file=fn,exist=iffile)
  if (iffile.eq. .true.) then
    isreal=.FALSE.;dmfile=fn;proceed=.TRUE.
  endif

  ! look for real density matrix file...
  inquire(file="density_matrix.inp",exist=iffile)
  if (iffile.eq. .true.) then
    isreal=.TRUE.;dmfile="density_matrix.inp";proceed=.TRUE.
  endif

  ! look for complex density matrix file...
  inquire(file="density_matrix_comp.inp",exist=iffile)
  if (iffile.eq. .true.) then
    isreal=.FALSE.;dmfile="density_matrix_comp.inp";proceed=.TRUE.
  endif


  if (proceed.eq. .true.) then
    print*,'***Atom', IATOM, '***'
    print*,"Reading from file ",dmfile
    open(65,file=dmfile)
    do ISP=1,NCDIJ
      tempdm=0.d0
      do l1=(ISP-1)*(2*lcor+1)+1,(2*lcor+1)*ISP
        ! different read statements for real vs. complex...
        if (isreal .eq. .TRUE.)then
          read(65,*)tempdm(l1-(ISP-1)*(2*lcor+1),:)
        else
          read(65,*)tempdm_c(l1-(ISP-1)*(2*lcor+1),:)
        endif
      enddo
      ! now we can assign the new density matrix that was read...
      if (isreal .eq. .TRUE.)then
        OCC_MAT(lcor**2+1:,lcor**2+1:,ISP)= tempdm(:,1:2*lcor+1)
      ! if complex then we need to add in the imaginary part...
      else
      OCC_MAT(lcor**2+1:,lcor**2+1:,ISP)= tempdm_c(:,1:2*lcor+1)
      OCC_MAT(lcor**2+1:,lcor**2+1:,ISP)=  &
          OCC_MAT(lcor**2+1:,lcor**2+1:,ISP) &
           +tempdm_c(:,2*lcor+2:4*lcor+2)*(0.0,1.0)
      endif
    enddo
    close(65)
  else
    print*,"Could not find a valid density matrix file to read"
  endif
endif

! END-MODVASP:1