src/gndstate.f90

   1
   2 ! Copyright (C) 2002-2007 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: gndstate
   8 ! !INTERFACE:
   9 subroutine gndstate
  10 ! !USES:
  11 use modmain
  12 ! !DESCRIPTION:
  13 !   Computes the self-consistent Kohn-Sham ground-state. General information is
  14 !   written to the file {\tt INFO.OUT}. First- and second-variational
  15 !   eigenvalues, eigenvectors and occupancies are written to the unformatted
  16 !   files {\tt EVALFV.OUT}, {\tt EVALSV.OUT}, {\tt EVECFV.OUT}, {\tt EVECSV.OUT}
  17 !   and {\tt OCCSV.OUT}.
  18 !
  19 ! !REVISION HISTORY:
  20 !   Created October 2002 (JKD)
  21 !EOP
  22 !BOC
  23 implicit none
  24 ! local variables
  25 logical exist
  26 integer ik,is,ia,idm
  27 integer n,nwork
  28 real(8) dv,timetot
  29 ! allocatable arrays
  30 real(8), allocatable :: v(:)
  31 real(8), allocatable :: work(:)
  32 real(8), allocatable :: evalfv(:,:)
  33 complex(8), allocatable :: evecfv(:,:,:)
  34 complex(8), allocatable :: evecsv(:,:)
  35 ! require forces for structural optimisation
  36 if ((task.eq.2).or.(task.eq.3)) tforce=.true.
  37 ! initialise global variables
  38 call init0
  39 call init1
  40 ! initialise OEP variables if required
  41 if (xctype.lt.0) call init2
  42 ! write the real and reciprocal lattice vectors to file
  43 call writelat
  44 ! write interatomic distances to file
  45 call writeiad(.false.)
  46 ! write symmetry matrices to file
  47 call writesym
  48 ! output the k-point set to file
  49 call writekpts
  50 ! write lattice vectors and atomic positions to file
  51 call writegeom(.false.)
  52 ! open INFO.OUT file
  53 open(60,file='INFO'//trim(filext),action='WRITE',form='FORMATTED')
  54 ! open TOTENERGY.OUT
  55 open(61,file='TOTENERGY'//trim(filext),action='WRITE',form='FORMATTED')
  56 ! open FERMIDOS.OUT
  57 open(62,file='FERMIDOS'//trim(filext),action='WRITE',form='FORMATTED')
  58 ! open MOMENT.OUT if required
  59 if (spinpol) open(63,file='MOMENT'//trim(filext),action='WRITE', &
  60  form='FORMATTED')
  61 ! open FORCEMAX.OUT if required
  62 if (tforce) open(64,file='FORCEMAX'//trim(filext),action='WRITE', &
  63  form='FORMATTED')
  64 ! open RMSDVEFF.OUT
  65 open(65,file='RMSDVEFF'//trim(filext),action='WRITE',form='FORMATTED')
  66 ! write out general information to INFO.OUT
  67 call writeinfo(60)
  68 ! initialise or read the charge density and potentials from file
  69 iscl=0
  70 write(60,*)
  71 if ((task.eq.1).or.(task.eq.3)) then
  72   call readstate
  73   write(60,'("Potential read in from STATE.OUT")')
  74 else if (task.eq.200) then
  75   call phveff
  76   write(60,'("Supercell potential constructed from STATE.OUT")')
  77 else
  78   call rhoinit
  79   call poteff
  80   call genveffig
  81   write(60,'("Density and potential initialised from atomic data")')
  82 end if
  83 call flushifc(60)
  84 ! size of mixing vector
  85 n=lmmaxvr*nrmtmax*natmtot+ngrtot
  86 if (spinpol) n=n*(1+ndmag)
  87 if (ldapu.ne.0) n=n+2*lmmaxlu*lmmaxlu*nspinor*nspinor*natmtot
  88 ! allocate mixing arrays
  89 allocate(v(n))
  90 ! determine the size of the mixer work array
  91 nwork=-1
  92 call mixerifc(mixtype,n,v,dv,nwork,v)
  93 allocate(work(nwork))
  94 ! set stop flag
  95 tstop=.false.
  96 10 continue
  97 ! set last iteration flag
  98 tlast=.false.
  99 ! delete any existing eigenvector files
 100 if ((task.eq.0).or.(task.eq.2)) call delevec
 101 ! begin the self-consistent loop
 102 write(60,*)
 103 write(60,'("+------------------------------+")')
 104 write(60,'("| Self-consistent loop started |")')
 105 write(60,'("+------------------------------+")')
 106 do iscl=1,maxscl
 107   write(60,*)
 108   write(60,'("+-------------------------+")')
 109   write(60,'("| Iteration number : ",I4," |")') iscl
 110   write(60,'("+-------------------------+")')
 111   call flushifc(60)
 112   if (iscl.ge.maxscl) then
 113     write(60,*)
 114     write(60,'("Reached self-consistent loops maximum")')
 115     tlast=.true.
 116   end if
 117   call flushifc(60)
 118 ! generate the core wavefunctions and densities
 119   call gencore
 120 ! find the new linearisation energies
 121   call linengy
 122 ! write out the linearisation energies
 123   call writelinen
 124 ! generate the APW radial functions
 125   call genapwfr
 126 ! generate the local-orbital radial functions
 127   call genlofr
 128 ! compute the overlap radial integrals
 129   call olprad
 130 ! compute the Hamiltonian radial integrals
 131   call hmlrad
 132 ! begin parallel loop over k-points
 133 !$OMP PARALLEL DEFAULT(SHARED) &
 134 !$OMP PRIVATE(evalfv,evecfv,evecsv)
 135 !$OMP DO
 136   do ik=1,nkpt
 137 ! every thread should allocate its own arrays
 138     allocate(evalfv(nstfv,nspnfv))
 139     allocate(evecfv(nmatmax,nstfv,nspnfv))
 140     allocate(evecsv(nstsv,nstsv))
 141 ! solve the first- and second-variational secular equations
 142     call seceqn(ik,evalfv,evecfv,evecsv)
 143 ! write the eigenvalues/vectors to file
 144     call putevalfv(ik,evalfv)
 145     call putevalsv(ik,evalsv(:,ik))
 146     call putevecfv(ik,evecfv)
 147     call putevecsv(ik,evecsv)
 148     deallocate(evalfv,evecfv,evecsv)
 149   end do
 150 !$OMP END DO
 151 !$OMP END PARALLEL
 152 ! find the occupation numbers and Fermi energy
 153   call occupy
 154 ! write out the eigenvalues and occupation numbers
 155   call writeeval
 156 ! write the Fermi energy to file
 157   call writefermi
 158 ! set the charge density and magnetisation to zero
 159   rhomt(:,:,:)=0.d0
 160   rhoir(:)=0.d0
 161   if (spinpol) then
 162     magmt(:,:,:,:)=0.d0
 163     magir(:,:)=0.d0
 164   end if
 165 !$OMP PARALLEL DEFAULT(SHARED) &
 166 !$OMP PRIVATE(evecfv,evecsv)
 167 !$OMP DO
 168   do ik=1,nkpt
 169     allocate(evecfv(nmatmax,nstfv,nspnfv))
 170     allocate(evecsv(nstsv,nstsv))
 171 ! write the occupancies to file
 172     call putoccsv(ik,occsv(:,ik))
 173 ! get the eigenvectors from file
 174     call getevecfv(vkl(:,ik),vgkl(:,:,:,ik),evecfv)
 175     call getevecsv(vkl(:,ik),evecsv)
 176 ! add to the density and magnetisation
 177     call rhovalk(ik,evecfv,evecsv)
 178     deallocate(evecfv,evecsv)
 179   end do
 180 !$OMP END DO
 181 !$OMP END PARALLEL
 182 ! symmetrise the density
 183   call symrf(lradstp,rhomt,rhoir)
 184 ! symmetrise the magnetisation
 185   if (spinpol) call symrvf(lradstp,magmt,magir)
 186 ! convert the density from a coarse to a fine radial mesh
 187   call rfmtctof(rhomt)
 188 ! convert the magnetisation from a coarse to a fine radial mesh
 189   do idm=1,ndmag
 190     call rfmtctof(magmt(:,:,:,idm))
 191   end do
 192 ! add the core density to the total density
 193   call addrhocr
 194 ! calculate the charges
 195   call charge
 196 ! calculate the moments
 197   if (spinpol) call moment
 198 ! normalise the density
 199   call rhonorm
 200 ! LDA+U
 201   if (ldapu.ne.0) then
 202 ! generate the LDA+U density matrix
 203     call gendmatlu
 204 ! generate the LDA+U potential matrix
 205     call genvmatlu
 206 ! write the LDA+U matrices to file
 207     call writeldapu
 208   end if
 209 ! compute the effective potential
 210   call poteff
 211 ! pack interstitial and muffin-tin effective potential and field into one array
 212   call packeff(.true.,n,v)
 213 ! mix in the old potential and field with the new
 214   call mixerifc(mixtype,n,v,dv,nwork,work)
 215 ! unpack potential and field
 216   call packeff(.false.,n,v)
 217 ! add the fixed spin moment effect field
 218   if (fixspin.ne.0) call fsmfield
 219 ! Fourier transform effective potential to G-space
 220   call genveffig
 221 ! reduce the external magnetic fields if required
 222   if (reducebf.lt.1.d0) then
 223     bfieldc(:)=bfieldc(:)*reducebf
 224     bfcmt(:,:,:)=bfcmt(:,:,:)*reducebf
 225   end if
 226 ! compute the energy components
 227   call energy
 228 ! output energy components
 229   call writeengy(60)
 230   write(60,*)
 231   write(60,'("Density of states at Fermi energy : ",G18.10)') fermidos
 232   write(60,'(" (states/Hartree/unit cell)")')
 233 ! write total energy to TOTENERGY.OUT and flush
 234   write(61,'(G22.12)') engytot
 235   call flushifc(61)
 236 ! write DOS at Fermi energy to FERMIDOS.OUT and flush
 237   write(62,'(G18.10)') fermidos
 238   call flushifc(62)
 239 ! output charges and moments
 240   call writechg(60)
 241 ! write total moment to MOMENT.OUT and flush
 242   if (spinpol) then
 243     write(63,'(3G18.10)') momtot(1:ndmag)
 244     call flushifc(63)
 245   end if
 246 ! output effective fields for fixed spin moment calculations
 247   if (fixspin.ne.0) call writefsm(60)
 248 ! check for WRITE file
 249   inquire(file='WRITE',exist=exist)
 250   if (exist) then
 251     write(60,*)
 252     write(60,'("WRITE file exists - writing STATE.OUT")')
 253     call writestate
 254     open(50,file='WRITE')
 255     close(50,status='DELETE')
 256   end if
 257 ! write STATE.OUT file if required
 258   if (nwrite.ge.1) then
 259     if (mod(iscl,nwrite).eq.0) then
 260       call writestate
 261       write(60,*)
 262       write(60,'("Wrote STATE.OUT")')
 263     end if
 264   end if
 265 ! exit self-consistent loop if last iteration is complete
 266   if (tlast) goto 20
 267 ! check for convergence
 268   if (iscl.ge.2) then
 269     write(60,*)
 270     write(60,'("RMS change in effective potential (target) : ",G18.10,&
 271      &" (",G18.10,")")') dv,epspot
 272     if (dv.lt.epspot) then
 273       write(60,*)
 274       write(60,'("Potential convergence target achieved")')
 275       tlast=.true.
 276     end if
 277     write(65,'(G18.10)') dv
 278     call flushifc(65)
 279   end if
 280   if (xctype.lt.0) then
 281     write(60,*)
 282     write(60,'("Magnitude of OEP residual : ",G18.10)') resoep
 283   end if
 284 ! check for STOP file
 285   inquire(file='STOP',exist=exist)
 286   if (exist) then
 287     write(60,*)
 288     write(60,'("STOP file exists - stopping self-consistent loop")')
 289     tstop=.true.
 290     tlast=.true.
 291     open(50,file='STOP')
 292     close(50,status='DELETE')
 293   end if
 294 ! output the current total CPU time
 295   timetot=timeinit+timemat+timefv+timesv+timerho+timepot+timefor
 296   write(60,*)
 297   write(60,'("Time (CPU seconds) : ",F12.2)') timetot
 298 ! end the self-consistent loop
 299 end do
 300 20 continue
 301 write(60,*)
 302 write(60,'("+------------------------------+")')
 303 write(60,'("| Self-consistent loop stopped |")')
 304 write(60,'("+------------------------------+")')
 305 ! write density and potentials to file only if maxscl > 1
 306 if (maxscl.gt.1) then
 307   call writestate
 308   write(60,*)
 309   write(60,'("Wrote STATE.OUT")')
 310 end if
 311 !-----------------------!
 312 !     compute forces    !
 313 !-----------------------!
 314 if ((.not.tstop).and.(tforce)) then
 315   call force
 316 ! output forces to INFO.OUT
 317   call writeforce(60)
 318 ! write maximum force magnitude to FORCEMAX.OUT
 319   write(64,'(G18.10)') forcemax
 320   call flushifc(64)
 321 end if
 322 !---------------------------------------!
 323 !     perform structural relaxation     !
 324 !---------------------------------------!
 325 if ((.not.tstop).and.((task.eq.2).or.(task.eq.3))) then
 326   write(60,*)
 327   write(60,'("Maximum force magnitude (target) : ",G18.10," (",G18.10,")")') &
 328    forcemax,epsforce
 329   call flushifc(60)
 330 ! check force convergence
 331   if (forcemax.le.epsforce) then
 332     write(60,*)
 333     write(60,'("Force convergence target achieved")')
 334     goto 30
 335   end if
 336 ! update the atomic positions if forces are not converged
 337   call updatpos
 338   write(60,*)
 339   write(60,'("+--------------------------+")')
 340   write(60,'("| Updated atomic positions |")')
 341   write(60,'("+--------------------------+")')
 342   do is=1,nspecies
 343     write(60,*)
 344     write(60,'("Species : ",I4," (",A,")")') is,trim(spsymb(is))
 345     write(60,'(" atomic positions (lattice) :")')
 346     do ia=1,natoms(is)
 347       write(60,'(I4," : ",3F14.8)') ia,atposl(:,ia,is)
 348     end do
 349   end do
 350 ! add blank line to TOTENERGY.OUT, FERMIDOS.OUT, MOMENT.OUT and RMSDVEFF.OUT
 351   write(61,*)
 352   write(62,*)
 353   if (spinpol) write (63,*)
 354   write(65,*)
 355 ! begin new self-consistent loop with updated positions
 356   goto 10
 357 end if
 358 30 continue
 359 ! output timing information
 360 write(60,*)
 361 write(60,'("Timings (CPU seconds) :")')
 362 write(60,'(" initialisation",T40,": ",F12.2)') timeinit
 363 write(60,'(" Hamiltonian and overlap matrix set up",T40,": ",F12.2)') timemat
 364 write(60,'(" first-variational secular equation",T40,": ",F12.2)') timefv
 365 if (spinpol) then
 366   write(60,'(" second-variational calculation",T40,": ",F12.2)') timesv
 367 end if
 368 write(60,'(" charge density calculation",T40,": ",F12.2)') timerho
 369 write(60,'(" potential calculation",T40,": ",F12.2)') timepot
 370 if (tforce) then
 371   write(60,'(" force calculation",T40,": ",F12.2)') timefor
 372 end if
 373 timetot=timeinit+timemat+timefv+timesv+timerho+timepot+timefor
 374 write(60,'(" total",T40,": ",F12.2)') timetot
 375 write(60,*)
 376 write(60,'("+----------------------------------+")')
 377 write(60,'("| EXCITING version ",I1.1,".",I1.1,".",I3.3," stopped |")') version
 378 write(60,'("+----------------------------------+")')
 379 ! close the INFO.OUT file
 380 close(60)
 381 ! close the TOTENERGY.OUT file
 382 close(61)
 383 ! close the FERMIDOS.OUT file
 384 close(62)
 385 ! close the MOMENT.OUT file
 386 if (spinpol) close(63)
 387 ! close the FORCEMAX.OUT file
 388 if (tforce) close(64)
 389 ! close the RMSDVEFF.OUT file
 390 close(65)
 391 deallocate(v,work)
 392 return
 393 end subroutine
 394 !EOC