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