src/init0.f90

   1
   2 ! Copyright (C) 2002-2005 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
   3 ! This file is distributed under the terms of the GNU General Public License.
   4 ! See the file COPYING for license details.
   5
   6 !BOP
   7 ! !ROUTINE: init0
   8 ! !INTERFACE:
   9 subroutine init0
  10 ! !USES:
  11 use modmain
  12 use modxcifc
  13 ! !DESCRIPTION:
  14 !   Performs basic consistency checks as well as allocating and initialising
  15 !   global variables not dependent on the $k$-point set.
  16 !
  17 ! !REVISION HISTORY:
  18 !   Created January 2004 (JKD)
  19 !EOP
  20 !BOC
  21 implicit none
  22 ! local variables
  23 integer is,js,ia,ias
  24 integer ist,l,m,lm,iv(3)
  25 real(8) ts0,ts1
  26
  27 !-------------------------------!
  28 !     zero timing variables     !
  29 !-------------------------------!
  30 timeinit=0.d0
  31 timemat=0.d0
  32 timefv=0.d0
  33 timesv=0.d0
  34 timerho=0.d0
  35 timepot=0.d0
  36 timefor=0.d0
  37 call timesec(ts0)
  38
  39 !------------------------------------!
  40 !     angular momentum variables     !
  41 !------------------------------------!
  42 lmmaxvr=(lmaxvr+1)**2
  43 lmmaxapw=(lmaxapw+1)**2
  44 lmmaxmat=(lmaxmat+1)**2
  45 lmmaxinr=(lmaxinr+1)**2
  46 if (lmaxvr.gt.lmaxapw) then
  47   write(*,*)
  48   write(*,'("Error(init0): lmaxvr > lmaxapw : ",2I8)') lmaxvr,lmaxapw
  49   write(*,*)
  50   stop
  51 end if
  52 if (lmaxmat.gt.lmaxapw) then
  53   write(*,*)
  54   write(*,'("Error(init0): lmaxmat > lmaxapw : ",2I8)') lmaxmat,lmaxapw
  55   write(*,*)
  56   stop
  57 end if
  58 ! index to (l,m) pairs
  59 if (allocated(idxlm)) deallocate(idxlm)
  60 allocate(idxlm(0:lmaxapw,-lmaxapw:lmaxapw))
  61 lm=0
  62 do l=0,lmaxapw
  63   do m=-l,l
  64     lm=lm+1
  65     idxlm(l,m)=lm
  66   end do
  67 end do
  68 ! array of i**l values
  69 if (allocated(zil)) deallocate(zil)
  70 allocate(zil(0:lmaxapw))
  71 do l=0,lmaxapw
  72   zil(l)=zi**l
  73 end do
  74
  75 !------------------------------------!
  76 !     index to atoms and species     !
  77 !------------------------------------!
  78 ! check if the system is an isolated molecule
  79 if (molecule) then
  80   primcell=.false.
  81   tshift=.false.
  82 end if
  83 ! find primitive cell if required
  84 if (primcell) call findprim
  85 natmmax=0
  86 ias=0
  87 do is=1,nspecies
  88   do ia=1,natoms(is)
  89     ias=ias+1
  90     idxas(ia,is)=ias
  91   end do
  92 ! maximum number of atoms over all species
  93   natmmax=max(natmmax,natoms(is))
  94 end do
  95 ! total number of atoms
  96 natmtot=ias
  97
  98 !------------------------!
  99 !     spin variables     !
 100 !------------------------!
 101 if (spinsprl) then
 102   select case(task)
 103   case(2,3,15,51,52,53,61,62,63,120,121)
 104     write(*,*)
 105     write(*,'("Error(init0): spin-spirals do not work with task ",I4)') task
 106     write(*,*)
 107     stop
 108   end select
 109   if (xctype.lt.0) then
 110     write(*,*)
 111     write(*,'("Error(init0): spin-spirals do not work with the OEP method")')
 112     write(*,*)
 113     stop
 114   end if
 115 end if
 116 ! spin-orbit coupling or fixed spin moment implies spin-polarised calculation
 117 if ((spinorb).or.(fixspin.ne.0).or.(spinsprl)) spinpol=.true.
 118 ! number of spinor components and maximum allowed occupancy
 119 if (spinpol) then
 120   nspinor=2
 121   occmax=1.d0
 122 else
 123   nspinor=1
 124   occmax=2.d0
 125 end if
 126 ! number of spin-dependent first-variational functions per state
 127 if (spinsprl) then
 128   nspnfv=2
 129 else
 130   nspnfv=1
 131 end if
 132 ! spin-polarised calculations require second-variational eigenvectors
 133 if (spinpol) tevecsv=.true.
 134 ! Hartree-Fock/RDMFT requires second-variational eigenvectors
 135 if ((task.eq.5).or.(task.eq.6).or.(task.eq.300)) tevecsv=.true.
 136 ! get exchange-correlation functional data
 137 call getxcdata(xctype,xcdescr,xcspin,xcgrad)
 138 if ((spinpol).and.(xcspin.eq.0)) then
 139   write(*,*)
 140   write(*,'("Error(init0): requested spin-polarised run with &
 141    &spin-unpolarised")')
 142   write(*,'(" exchange-correlation functional")')
 143   write(*,*)
 144   stop
 145 end if
 146 ! check for collinearity in the z-direction and set the dimension of the
 147 ! magnetisation and exchange-correlation vector fields
 148 if (spinpol) then
 149   ndmag=1
 150   if ((abs(bfieldc(1)).gt.epslat).or.(abs(bfieldc(2)).gt.epslat)) ndmag=3
 151   do is=1,nspecies
 152     do ia=1,natoms(is)
 153       if ((abs(bfcmt(1,ia,is)).gt.epslat).or.(abs(bfcmt(2,ia,is)).gt.epslat)) &
 154        ndmag=3
 155     end do
 156   end do
 157 ! source-free fields and spin-spirals are non-collinear in general
 158   if ((nosource).or.(spinsprl)) ndmag=3
 159 ! spin-orbit coupling is non-collinear in general
 160   if (spinorb) ndmag=3
 161 else
 162   ndmag=0
 163 end if
 164 ! set the non-collinear flag
 165 if (ndmag.eq.3) then
 166   ncmag=.true.
 167 else
 168   ncmag=.false.
 169 end if
 170 if ((ncmag).and.(xcgrad.gt.0)) then
 171   write(*,*)
 172   write(*,'("Warning(init0): GGA inconsistent with non-collinear magnetism")')
 173 end if
 174 ! set fixed spin moment effective field to zero
 175 bfsmc(:)=0.d0
 176 ! set muffin-tin FSM fields to zero
 177 bfsmcmt(:,:,:)=0.d0
 178
 179 !-------------------------------------!
 180 !     lattice and symmetry set up     !
 181 !-------------------------------------!
 182 ! generate the reciprocal lattice vectors and unit cell volume
 183 call reciplat
 184 ! compute the inverse of the lattice vector matrix
 185 call r3minv(avec,ainv)
 186 ! compute the inverse of the reciprocal vector matrix
 187 call r3minv(bvec,binv)
 188 do is=1,nspecies
 189   do ia=1,natoms(is)
 190 ! map atomic lattice coordinates to [0,1) if not in molecule mode
 191     if (.not.molecule) call r3frac(epslat,atposl(:,ia,is),iv)
 192 ! determine atomic Cartesian coordinates
 193     call r3mv(avec,atposl(:,ia,is),atposc(:,ia,is))
 194 ! lattice coordinates of the muffin-tin magnetic fields
 195     call r3mv(ainv,bfcmt(:,ia,is),bflmt(:,ia,is))
 196   end do
 197 end do
 198 ! lattice coordinates of the global magnetic field
 199 call r3mv(ainv,bfieldc,bfieldl)
 200 ! Cartesian coordinates of the spin-spiral vector
 201 call r3mv(bvec,vqlss,vqcss)
 202 ! find Bravais lattice symmetries
 203 call findsymlat
 204 ! use only the identity if required
 205 if (nosym) nsymlat=1
 206 ! find the crystal symmetries and shift atomic positions if required
 207 call findsymcrys
 208 ! find the site symmetries
 209 call findsymsite
 210 ! automatically determine the muffin-tin radii if required
 211 if (autormt) call autoradmt
 212 ! check for overlapping muffin-tins
 213 call checkmt
 214
 215 !-----------------------!
 216 !     radial meshes     !
 217 !-----------------------!
 218 nrmtmax=1
 219 nrcmtmax=1
 220 js=1
 221 do is=1,nspecies
 222 ! make the muffin-tin mesh commensurate with lradstp
 223   nrmt(is)=nrmt(is)-mod(nrmt(is)-1,lradstp)
 224   nrmtmax=max(nrmtmax,nrmt(is))
 225 ! number of coarse radial mesh points
 226   nrcmt(is)=(nrmt(is)-1)/lradstp+1
 227   nrcmtmax=max(nrcmtmax,nrcmt(is))
 228 ! smallest muffin-tin radius
 229   if (rmt(is).lt.rmt(js)) js=is
 230 end do
 231 if ((isgkmax.lt.1).or.(isgkmax.gt.nspecies)) isgkmax=js
 232 ! set up atomic and muffin-tin radial meshes
 233 call genrmesh
 234
 235 !--------------------------------------!
 236 !     charges and number of states     !
 237 !--------------------------------------!
 238 chgzn=0.d0
 239 chgcr=0.d0
 240 chgval=0.d0
 241 spnstmax=0
 242 do is=1,nspecies
 243 ! nuclear charge
 244   chgzn=chgzn+spzn(is)*dble(natoms(is))
 245 ! find the maximum number of atomic states
 246   spnstmax=max(spnstmax,spnst(is))
 247 ! compute the electronic charge for each species, as well as the total core and
 248 ! valence charge
 249   spze(is)=0.d0
 250   do ist=1,spnst(is)
 251     spze(is)=spze(is)+spocc(ist,is)
 252     if (spcore(ist,is)) then
 253       chgcr=chgcr+dble(natoms(is))*spocc(ist,is)
 254     else
 255       chgval=chgval+dble(natoms(is))*spocc(ist,is)
 256     end if
 257   end do
 258 end do
 259 ! add excess charge
 260 chgval=chgval+chgexs
 261 ! total charge
 262 chgtot=chgcr+chgval
 263 if (chgtot.lt.1.d-8) then
 264   write(*,*)
 265   write(*,'("Error(init0): zero total charge")')
 266   write(*,*)
 267   stop
 268 end if
 269 ! effective Wigner radius
 270 rwigner=(3.d0/(fourpi*(chgtot/omega)))**(1.d0/3.d0)
 271
 272 !-------------------------!
 273 !     G-vector arrays     !
 274 !-------------------------!
 275 ! find the G-vector grid sizes
 276 call gridsize
 277 ! generate the G-vectors
 278 call gengvec
 279 ! generate the spherical harmonics of the G-vectors
 280 call genylmg
 281 ! allocate structure factor array for G-vectors
 282 if (allocated(sfacg)) deallocate(sfacg)
 283 allocate(sfacg(ngvec,natmtot))
 284 ! generate structure factors for G-vectors
 285 call gensfacgp(ngvec,vgc,ngvec,sfacg)
 286 ! generate the characteristic function
 287 call gencfun
 288
 289 !-------------------------!
 290 !     atoms and cores     !
 291 !-------------------------!
 292 ! solve the Kohn-Sham-Dirac equations for all atoms
 293 call allatoms
 294 ! allocate core state eigenvalue array and set to default
 295 if (allocated(evalcr)) deallocate(evalcr)
 296 allocate(evalcr(spnstmax,natmtot))
 297 do is=1,nspecies
 298   do ia=1,natoms(is)
 299     ias=idxas(ia,is)
 300     do ist=1,spnst(is)
 301       evalcr(ist,ias)=speval(ist,is)
 302     end do
 303   end do
 304 end do
 305 ! allocate core state radial wavefunction array
 306 if (allocated(rwfcr)) deallocate(rwfcr)
 307 allocate(rwfcr(spnrmax,2,spnstmax,natmtot))
 308 ! allocate core state charge density array
 309 if (allocated(rhocr)) deallocate(rhocr)
 310 allocate(rhocr(spnrmax,natmtot))
 311
 312 !---------------------------------------!
 313 !     charge density and potentials     !
 314 !---------------------------------------!
 315 ! allocate charge density arrays
 316 if (allocated(rhomt)) deallocate(rhomt)
 317 allocate(rhomt(lmmaxvr,nrmtmax,natmtot))
 318 if (allocated(rhoir)) deallocate(rhoir)
 319 allocate(rhoir(ngrtot))
 320 ! allocate magnetisation arrays
 321 if (allocated(magmt)) deallocate(magmt)
 322 if (allocated(magir)) deallocate(magir)
 323 if (spinpol) then
 324   allocate(magmt(lmmaxvr,nrmtmax,natmtot,ndmag))
 325   allocate(magir(ngrtot,ndmag))
 326 end if
 327 ! Coulomb potential
 328 if (allocated(vclmt)) deallocate(vclmt)
 329 allocate(vclmt(lmmaxvr,nrmtmax,natmtot))
 330 if (allocated(vclir)) deallocate(vclir)
 331 allocate(vclir(ngrtot))
 332 ! exchange-correlation potential
 333 if (allocated(vxcmt)) deallocate(vxcmt)
 334 allocate(vxcmt(lmmaxvr,nrmtmax,natmtot))
 335 if (allocated(vxcir)) deallocate(vxcir)
 336 allocate(vxcir(ngrtot))
 337 ! exchange-correlation magnetic field
 338 if (allocated(bxcmt)) deallocate(bxcmt)
 339 if (allocated(bxcir)) deallocate(bxcir)
 340 if (spinpol) then
 341   allocate(bxcmt(lmmaxvr,nrmtmax,natmtot,ndmag))
 342   allocate(bxcir(ngrtot,ndmag))
 343 end if
 344 ! exchange energy density
 345 if (allocated(exmt)) deallocate(exmt)
 346 allocate(exmt(lmmaxvr,nrmtmax,natmtot))
 347 if (allocated(exir)) deallocate(exir)
 348 allocate(exir(ngrtot))
 349 ! correlation energy density
 350 if (allocated(ecmt)) deallocate(ecmt)
 351 allocate(ecmt(lmmaxvr,nrmtmax,natmtot))
 352 if (allocated(ecir)) deallocate(ecir)
 353 allocate(ecir(ngrtot))
 354 ! effective potential
 355 if (allocated(veffmt)) deallocate(veffmt)
 356 allocate(veffmt(lmmaxvr,nrmtmax,natmtot))
 357 if (allocated(veffir)) deallocate(veffir)
 358 allocate(veffir(ngrtot))
 359 if (allocated(veffig)) deallocate(veffig)
 360 allocate(veffig(ngvec))
 361 ! allocate muffin-tin charge and moment arrays
 362 if (allocated(chgmt)) deallocate(chgmt)
 363 allocate(chgmt(natmtot))
 364 if (allocated(mommt)) deallocate(mommt)
 365 allocate(mommt(3,natmtot))
 366
 367 !--------------------------------------------!
 368 !     forces and structural optimisation     !
 369 !--------------------------------------------!
 370 if (allocated(forcehf)) deallocate(forcehf)
 371 allocate(forcehf(3,natmtot))
 372 if (allocated(forcecr)) deallocate(forcecr)
 373 allocate(forcecr(3,natmtot))
 374 if (allocated(forceibs)) deallocate(forceibs)
 375 allocate(forceibs(3,natmtot))
 376 if (allocated(forcetot)) deallocate(forcetot)
 377 allocate(forcetot(3,natmtot))
 378 if (allocated(forcetp)) deallocate(forcetp)
 379 allocate(forcetp(3,natmtot))
 380 if (allocated(tauatm)) deallocate(tauatm)
 381 allocate(tauatm(natmtot))
 382 ! initialise the previous force
 383 forcetp(:,:)=0.d0
 384 ! initial step sizes
 385 tauatm(:)=tau0atm
 386
 387 !-------------------------!
 388 !     LDA+U variables     !
 389 !-------------------------!
 390 if ((ldapu.ne.0).or.(task.eq.17)) then
 391 ! LDA+U requires second-variational eigenvectors
 392   tevecsv=.true.
 393 ! density matrices
 394   if (allocated(dmatlu)) deallocate(dmatlu)
 395   allocate(dmatlu(lmmaxlu,lmmaxlu,nspinor,nspinor,natmtot))
 396 ! potential matrix elements
 397   if (allocated(vmatlu)) deallocate(vmatlu)
 398   allocate(vmatlu(lmmaxlu,lmmaxlu,nspinor,nspinor,natmtot))
 399 ! zero the potential
 400   vmatlu(:,:,:,:,:)=0.d0
 401 ! energy for each atom
 402   if (allocated(engyalu)) deallocate(engyalu)
 403   allocate(engyalu(natmtot))
 404 ! interpolation constants (alpha)
 405   if (allocated(alphalu)) deallocate(alphalu)
 406   allocate(alphalu(natmtot))
 407 end if
 408
 409 !-----------------------!
 410 !     miscellaneous     !
 411 !-----------------------!
 412 ! determine the nuclear-nuclear energy
 413 call energynn
 414 ! get smearing function data
 415 call getsdata(stype,sdescr)
 416 ! generate the spherical harmonic transform (SHT) matrices
 417 call genshtmat
 418 ! allocate 1D plotting arrays
 419 if (allocated(dvp1d)) deallocate(dvp1d)
 420 allocate(dvp1d(nvp1d))
 421 if (allocated(vplp1d)) deallocate(vplp1d)
 422 allocate(vplp1d(3,npp1d))
 423 if (allocated(dpp1d)) deallocate(dpp1d)
 424 allocate(dpp1d(npp1d))
 425 ! zero self-consistent loop number
 426 iscl=0
 427 tlast=.false.
 428
 429 call timesec(ts1)
 430 timeinit=timeinit+ts1-ts0
 431
 432 return
 433 end subroutine
 434 !EOC
 435