src/readinput.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: readinput
   8 ! !INTERFACE:
   9 subroutine readinput
  10 ! !USES:
  11 use modmain
  12 ! !DESCRIPTION:
  13 !   Reads in the input parameters from the file {\tt exciting.in} as well as
  14 !   from the species files. Also sets default values for the input parameters.
  15 !
  16 ! !REVISION HISTORY:
  17 !   Created September 2002 (JKD)
  18 !EOP
  19 !BOC
  20 implicit none
  21 ! local variables
  22 integer is,js,ia,ja,ias
  23 integer i,l,iv,iostat
  24 integer ist,io,nlx,ilx,lx,ilo
  25 real(8) sc,sc1,sc2,sc3
  26 real(8) vacuum,v(3),t1,t2
  27 character(256) fname
  28 character(256) str
  29 character(256) bname
  30 character(256) sppath
  31
  32 !------------------------!
  33 !     default values     !
  34 !------------------------!
  35 ntasks=1
  36 tasks(1)=-1
  37 avec(:,:)=0.d0
  38 avec(1,1)=1.d0
  39 avec(2,2)=1.d0
  40 avec(3,3)=1.d0
  41 sc=1.d0
  42 sc1=1.d0
  43 sc2=1.d0
  44 sc3=1.d0
  45 epslat=1.d-6
  46 primcell=.false.
  47 tshift=.true.
  48 ngridk(:)=1
  49 vkloff(:)=0.d0
  50 rlambda=30.d0
  51 autokpt=.false.
  52 reducek=.true.
  53 ngridq(:)=1
  54 reduceq=.true.
  55 rgkmax=7.d0
  56 gmaxvr=12.d0
  57 lmaxapw=8
  58 lmaxvr=7
  59 lmaxmat=5
  60 lmaxinr=2
  61 fracinr=0.25d0
  62 npsden=9
  63 xctype=3
  64 stype=0
  65 swidth=0.01d0
  66 epsocc=1.d-8
  67 epschg=1.d-3
  68 nempty=5
  69 beta0=0.1d0
  70 betamax=1.d0
  71 maxscl=200
  72 epspot=1.d-6
  73 epsengy=1.d-7
  74 epsforce=5.d-4
  75 cfdamp=0.d0
  76 molecule=.false.
  77 vacuum=10.d0
  78 nspecies=0
  79 natoms(:)=0
  80 sppath='./'
  81 scrpath='./'
  82 nvp1d=2
  83 if (allocated(vvlp1d)) deallocate(vvlp1d)
  84 allocate(vvlp1d(3,nvp1d))
  85 vvlp1d(:,1)=0.d0
  86 vvlp1d(:,2)=1.d0
  87 npp1d=200
  88 vclp2d(:,:)=0.d0
  89 vclp2d(1,2)=1.d0
  90 vclp2d(2,3)=1.d0
  91 np2d(:)=40
  92 nup3d(:)=1
  93 np3d(:)=20
  94 nwdos=500
  95 ngrdos=100
  96 nsmdos=0
  97 wdos(1)=-0.5d0
  98 wdos(2)=0.5d0
  99 bcsym=.true.
 100 spinpol=.false.
 101 spinorb=.false.
 102 tau0atm=0.2d0
 103 nstfsp=6
 104 lradstp=4
 105 chgexs=0.d0
 106 nprad=4
 107 scissor=0.d0
 108 noptcomp=1
 109 optcomp(:,1)=1
 110 usegdft=.false.
 111 intraband=.false.
 112 evalmin=-4.5d0
 113 deband=0.0025d0
 114 bfieldc(:)=0.d0
 115 fixspin=0
 116 momfix(:)=0.d0
 117 mommtfix(:,:,:)=0.d0
 118 taufsm=0.01d0
 119 autormt=.false.
 120 rmtapm(1)=0.25d0
 121 rmtapm(2)=0.95d0
 122 nosym=.false.
 123 deltaph=0.03d0
 124 nphwrt=1
 125 if (allocated(vqlwrt)) deallocate(vqlwrt)
 126 allocate(vqlwrt(3,nphwrt))
 127 vqlwrt(:,:)=0.d0
 128 notelns=0
 129 tforce=.false.
 130 tfibs=.true.
 131 maxitoep=120
 132 tauoep(1)=1.d0
 133 tauoep(2)=0.2d0
 134 tauoep(3)=1.5d0
 135 nkstlist=1
 136 kstlist(:,1)=1
 137 vklem(:)=0.d0
 138 deltaem=0.025d0
 139 ndspem=1
 140 nosource=.false.
 141 spinsprl=.false.
 142 vqlss(:)=0.d0
 143 nwrite=0
 144 tevecsv=.false.
 145 ldapu=0
 146 llu(:)=-1
 147 ujlu(:,:)=0.d0
 148 rdmxctype=2
 149 rdmmaxscl=1
 150 maxitn=250
 151 maxitc=10
 152 taurdmn=1.d0
 153 taurdmc=0.5d0
 154 rdmalpha=0.7d0
 155 reducebf=1.d0
 156
 157 !-------------------------------!
 158 !     read from exciting.in     !
 159 !-------------------------------!
 160 fname='exciting.in'
 161 !----- uncomment for command line input filename (M. Rajagopalan) --------------
 162 !if (iargc()>1) then
 163 !  write(*,*)
 164 !  write(*,'("Usage: exciting [INPUT FILE]")')
 165 !  write(*,*)
 166 !  stop
 167 !end if
 168 !if (iargc().eq.1) then
 169 !  call getarg(1,fname)
 170 !end if
 171 !-------------------------------------------------------------------------------
 172 open(50,file=trim(fname),action='READ',status='OLD',form='FORMATTED', &
 173  iostat=iostat)
 174 if (iostat.ne.0) then
 175   write(*,*)
 176   write(*,'("Error(readinput): error opening ",A)') trim(fname)
 177   write(*,*)
 178   stop
 179 end if
 180 10 continue
 181 read(50,*,end=30) bname
 182 ! check for a comment
 183 if ((scan(trim(bname),'!').eq.1).or.(scan(trim(bname),'#').eq.1)) goto 10
 184 select case(trim(bname))
 185 case('tasks')
 186   do i=1,maxtasks
 187     read(50,'(A80)') str
 188     if (trim(str).eq.'') then
 189       if (i.eq.1) then
 190         write(*,*)
 191         write(*,'("Error(readinput): no tasks to perform")')
 192         write(*,*)
 193         stop
 194       end if
 195       ntasks=i-1
 196       goto 10
 197     end if
 198     read(str,*,iostat=iostat) tasks(i)
 199     if (iostat.ne.0) then
 200       write(*,*)
 201       write(*,'("Error(readinput): error reading tasks")')
 202       write(*,'("(blank line required after tasks block)")')
 203       write(*,*)
 204       stop
 205     end if
 206   end do
 207   write(*,*)
 208   write(*,'("Error(readinput): too many tasks")')
 209   write(*,*)
 210   stop
 211 case('avec')
 212   read(50,*,err=20) avec(1,1),avec(2,1),avec(3,1)
 213   read(50,*,err=20) avec(1,2),avec(2,2),avec(3,2)
 214   read(50,*,err=20) avec(1,3),avec(2,3),avec(3,3)
 215 case('scale')
 216   read(50,*,err=20) sc
 217 case('scale1')
 218   read(50,*,err=20) sc1
 219 case('scale2')
 220   read(50,*,err=20) sc2
 221 case('scale3')
 222   read(50,*,err=20) sc3
 223 case('epslat')
 224   read(50,*,err=20) epslat
 225   if (epslat.le.0.d0) then
 226     write(*,*)
 227     write(*,'("Error(readinput): epslat <= 0 : ",G18.10)') epslat
 228     write(*,*)
 229     stop
 230   end if
 231 case('primcell')
 232   read(50,*,err=20) primcell
 233 case('tshift')
 234   read(50,*,err=20) tshift
 235 case('rlambda')
 236   read(50,*,err=20) rlambda
 237   if (rlambda.le.0.d0) then
 238     write(*,*)
 239     write(*,'("Error(readinput): rlambda <= 0 : ",G18.10)') rlambda
 240     write(*,*)
 241     stop
 242   end if
 243 case('autokpt')
 244   read(50,*,err=20) autokpt
 245 case('ngridk')
 246   read(50,*,err=20) ngridk(1),ngridk(2),ngridk(3)
 247   if ((ngridk(1).le.0).or.(ngridk(2).le.0).or.(ngridk(3).le.0)) then
 248     write(*,*)
 249     write(*,'("Error(readinput): invalid ngridk : ",3I8)') ngridk
 250     write(*,*)
 251     stop
 252   end if
 253 case('vkloff')
 254   read(50,*,err=20) vkloff(1),vkloff(2),vkloff(3)
 255 case('reducek')
 256   read(50,*,err=20) reducek
 257 case('ngridq')
 258   read(50,*,err=20) ngridq(1),ngridq(2),ngridq(3)
 259   if ((ngridq(1).le.0).or.(ngridq(2).le.0).or.(ngridq(3).le.0)) then
 260     write(*,*)
 261     write(*,'("Error(readinput): invalid ngridq : ",3I8)') ngridq
 262     write(*,*)
 263     stop
 264   end if
 265 case('reduceq')
 266   read(50,*,err=20) reduceq
 267 case('rgkmax')
 268   read(50,*,err=20) rgkmax
 269   if (rgkmax.le.0.d0) then
 270     write(*,*)
 271     write(*,'("Error(readinput): rgkmax <= 0 : ",G18.10)') rgkmax
 272     write(*,*)
 273     stop
 274   end if
 275 case('gmaxvr')
 276   read(50,*,err=20) gmaxvr
 277 case('lmaxapw')
 278   read(50,*,err=20) lmaxapw
 279   if (lmaxapw.lt.0) then
 280     write(*,*)
 281     write(*,'("Error(readinput): lmaxapw < 0 : ",I8)') lmaxapw
 282     write(*,*)
 283     stop
 284   end if
 285   if (lmaxapw.ge.maxlapw) then
 286     write(*,*)
 287     write(*,'("Error(readinput): lmaxapw too large : ",I8)') lmaxapw
 288     write(*,'("Adjust maxlapw in modmain and recompile code")')
 289     write(*,*)
 290     stop
 291   end if
 292 case('lmaxvr')
 293   read(50,*,err=20) lmaxvr
 294   if (lmaxvr.lt.3) then
 295     write(*,*)
 296     write(*,'("Error(readinput): lmaxvr < 3 : ",I8)') lmaxvr
 297     write(*,*)
 298     stop
 299   end if
 300 case('lmaxmat')
 301   read(50,*,err=20) lmaxmat
 302   if (lmaxmat.lt.0) then
 303     write(*,*)
 304     write(*,'("Error(readinput): lmaxmat < 0 : ",I8)') lmaxmat
 305     write(*,*)
 306     stop
 307   end if
 308 case('lmaxinr')
 309   read(50,*,err=20) lmaxinr
 310   if (lmaxinr.lt.0) then
 311     write(*,*)
 312     write(*,'("Error(readinput): lmaxinr < 0 : ",I8)') lmaxinr
 313     write(*,*)
 314     stop
 315   end if
 316 case('fracinr')
 317   read(50,*,err=20) fracinr
 318 case('npsden')
 319   read(50,*,err=20) npsden
 320   if (npsden.lt.2) then
 321     write(*,*)
 322     write(*,'("Error(readinput): npsden < 2 : ",I8)') npsden
 323     write(*,*)
 324     stop
 325   end if
 326 case('spinpol')
 327   read(50,*,err=20) spinpol
 328 case('spinorb')
 329   read(50,*,err=20) spinorb
 330 case('xctype')
 331   read(50,*,err=20) xctype
 332 case('stype')
 333   read(50,*,err=20) stype
 334 case('swidth')
 335   read(50,*,err=20) swidth
 336   if (swidth.lt.1.d-9) then
 337     write(*,*)
 338     write(*,'("Error(readinput): swidth too small or negative : ",G18.10)') &
 339      swidth
 340     write(*,*)
 341     stop
 342   end if
 343 case('epsocc')
 344   read(50,*,err=20) epsocc
 345   if (epsocc.le.0.d0) then
 346     write(*,*)
 347     write(*,'("Error(readinput): epsocc <= 0 : ",G18.10)') epsocc
 348     write(*,*)
 349     stop
 350   end if
 351 case('epschg')
 352   read(50,*,err=20) epschg
 353   if (epschg.le.0.d0) then
 354     write(*,*)
 355     write(*,'("Error(readinput): epschg <= 0 : ",G18.10)') epschg
 356     write(*,*)
 357     stop
 358   end if
 359 case('nempty')
 360   read(50,*,err=20) nempty
 361   if (nempty.le.0) then
 362     write(*,*)
 363     write(*,'("Error(readinput): nempty <= 0 : ",I8)') nempty
 364     write(*,*)
 365     stop
 366   end if
 367 case('beta0')
 368   read(50,*,err=20) beta0
 369   if (beta0.lt.0.d0) then
 370     write(*,*)
 371     write(*,'("Error(readinput): beta0 < 0 : ",G18.10)') beta0
 372     write(*,*)
 373     stop
 374   end if
 375 case('betamax')
 376   read(50,*,err=20) betamax
 377   if (betamax.lt.0.d0) then
 378     write(*,*)
 379     write(*,'("Error(readinput): betamax < 0 : ",G18.10)') betamax
 380     write(*,*)
 381     stop
 382   end if
 383 case('maxscl')
 384   read(50,*,err=20) maxscl
 385   if (maxscl.lt.0) then
 386     write(*,*)
 387     write(*,'("Error(readinput): maxscl < 0 : ",I8)') maxscl
 388     write(*,*)
 389     stop
 390   end if
 391 case('epspot')
 392   read(50,*,err=20) epspot
 393 case('epsengy')
 394   read(50,*,err=20) epsengy
 395 case('epsforce')
 396   read(50,*,err=20) epsforce
 397 case('cfdamp')
 398   read(50,*,err=20) cfdamp
 399   if (cfdamp.lt.0.d0) then
 400     write(*,*)
 401     write(*,'("Error(readinput): cfdamp < 0 : ",G18.10)') cfdamp
 402     write(*,*)
 403     stop
 404   end if
 405 case('sppath')
 406   read(50,*,err=20) sppath
 407   sppath=adjustl(sppath)
 408 case('scrpath')
 409   read(50,*,err=20) scrpath
 410 case('molecule')
 411   read(50,*,err=20) molecule
 412 case('vacuum')
 413   read(50,*,err=20) vacuum
 414   if (vacuum.lt.0.d0) then
 415     write(*,*)
 416     write(*,'("Error(readinput): vacuum < 0 : ",G18.10)') vacuum
 417     write(*,*)
 418     stop
 419   end if
 420 case('atoms')
 421   read(50,*,err=20) nspecies
 422   if (nspecies.le.0) then
 423     write(*,*)
 424     write(*,'("Error(readinput): nspecies <= 0 : ",I8)') nspecies
 425     write(*,*)
 426     stop
 427   end if
 428   if (nspecies.gt.maxspecies) then
 429     write(*,*)
 430     write(*,'("Error(readinput): nspecies too large : ",I8)') nspecies
 431     write(*,'("Adjust maxspecies in modmain and recompile code")')
 432     write(*,*)
 433     stop
 434   end if
 435   do is=1,nspecies
 436     read(50,*,err=20) spfname(is)
 437     spfname(is)=adjustl(spfname(is))
 438     read(50,*,err=20) natoms(is)
 439     if (natoms(is).le.0) then
 440       write(*,*)
 441       write(*,'("Error(readinput): natoms <= 0 : ",I8)') natoms(is)
 442       write(*,'(" for species ",I4)') is
 443       write(*,*)
 444       stop
 445     end if
 446     if (natoms(is).gt.maxatoms) then
 447       write(*,*)
 448       write(*,'("Error(readinput): natoms too large : ",I8)') natoms(is)
 449       write(*,'(" for species ",I4)') is
 450       write(*,'("Adjust maxatoms in modmain and recompile code")')
 451       write(*,*)
 452       stop
 453     end if
 454     do ia=1,natoms(is)
 455       read(50,*,err=20) atposl(1,ia,is),atposl(2,ia,is),atposl(3,ia,is), &
 456        bfcmt(1,ia,is),bfcmt(2,ia,is),bfcmt(3,ia,is)
 457     end do
 458   end do
 459 case('plot1d')
 460   read(50,*,err=20) nvp1d,npp1d
 461   if (nvp1d.lt.1) then
 462     write(*,*)
 463     write(*,'("Error(readinput): nvp1d < 1 : ",I8)') nvp1d
 464     write(*,*)
 465     stop
 466   end if
 467   if (npp1d.lt.nvp1d) then
 468     write(*,*)
 469     write(*,'("Error(readinput): npp1d < nvp1d : ",2I8)') npp1d,nvp1d
 470     write(*,*)
 471     stop
 472   end if
 473   if (allocated(vvlp1d)) deallocate(vvlp1d)
 474   allocate(vvlp1d(3,nvp1d))
 475   do iv=1,nvp1d
 476     read(50,*,err=20) vvlp1d(1,iv),vvlp1d(2,iv),vvlp1d(3,iv)
 477   end do
 478 case('plot2d')
 479   read(50,*,err=20) vclp2d(1,1),vclp2d(2,1),vclp2d(3,1)
 480   read(50,*,err=20) vclp2d(1,2),vclp2d(2,2),vclp2d(3,2)
 481   read(50,*,err=20) vclp2d(1,3),vclp2d(2,3),vclp2d(3,3)
 482   read(50,*,err=20) np2d(1),np2d(2)
 483   if ((np2d(1).lt.1).or.(np2d(2).lt.1)) then
 484     write(*,*)
 485     write(*,'("Error(readinput): np2d < 1 : ",2I8)') np2d
 486     write(*,*)
 487     stop
 488   end if
 489 case('plot3d')
 490   read(50,*,err=20) nup3d(1),nup3d(2),nup3d(3)
 491   if ((nup3d(1).lt.1).or.(nup3d(2).lt.1).or.(nup3d(3).lt.1)) then
 492     write(*,*)
 493     write(*,'("Error(readinput): nup3d < 1 : ",3I8)') nup3d
 494     write(*,*)
 495     stop
 496   end if
 497   read(50,*,err=20) np3d(1),np3d(2),np3d(3)
 498   if ((np3d(1).lt.1).or.(np3d(2).lt.1).or.(np3d(3).lt.1)) then
 499     write(*,*)
 500     write(*,'("Error(readinput): np3d < 1 : ",3I8)') np3d
 501     write(*,*)
 502     stop
 503   end if
 504 case('dos')
 505   read(50,*,err=20) nwdos,ngrdos,nsmdos
 506   if (nwdos.lt.2) then
 507     write(*,*)
 508     write(*,'("Error(readinput): nwdos < 2 : ",I8)') nwdos
 509     write(*,*)
 510     stop
 511   end if
 512   if (ngrdos.lt.1) then
 513     write(*,*)
 514     write(*,'("Error(readinput): ngrdos < 1 : ",I8)') ngrdos
 515     write(*,*)
 516     stop
 517   end if
 518   if (nsmdos.lt.0) then
 519     write(*,*)
 520     write(*,'("Error(readinput): nsmdos < 0 : ",I8)') nsmdos
 521     write(*,*)
 522     stop
 523   end if
 524   read(50,*,err=20) wdos(1),wdos(2)
 525   if (wdos(1).ge.wdos(2)) then
 526     write(*,*)
 527     write(*,'("Error(readinput): wdos(1) >= wdos(2) : ",2G18.10)') wdos
 528     write(*,*)
 529     stop
 530   end if
 531 case('bcsym')
 532   read(50,*,err=20) bcsym
 533 case('tau0atm')
 534   read(50,*,err=20) tau0atm
 535 case('nstfsp')
 536   read(50,*,err=20) nstfsp
 537   if (nstfsp.le.0) then
 538     write(*,*)
 539     write(*,'("Error(readinput): nstfsp <= 0 : ",I8)') nstfsp
 540     write(*,*)
 541     stop
 542   end if
 543 case('lradstp')
 544   read(50,*,err=20) lradstp
 545   if (lradstp.le.0) then
 546     write(*,*)
 547     write(*,'("Error(readinput): lradstp <= 0 : ",I8)') lradstp
 548     write(*,*)
 549     stop
 550   end if
 551 case('chgexs')
 552   read(50,*,err=20) chgexs
 553 case('nprad')
 554   read(50,*,err=20) nprad
 555   if (nprad.lt.2) then
 556     write(*,*)
 557     write(*,'("Error(readinput): nprad < 2 : ",I8)') nprad
 558     write(*,*)
 559     stop
 560   end if
 561 case('scissor')
 562   read(50,*,err=20) scissor
 563 case('optcomp')
 564   do i=1,27
 565     read(50,'(A80)') str
 566     if (trim(str).eq.'') then
 567       if (i.eq.1) then
 568         write(*,*)
 569         write(*,'("Error(readinput): empty optical component list")')
 570         write(*,*)
 571         stop
 572       end if
 573       noptcomp=i-1
 574       goto 10
 575     end if
 576     read(str,*,iostat=iostat) optcomp(:,i)
 577     if (iostat.ne.0) then
 578       write(*,*)
 579       write(*,'("Error(readinput): error reading optical component list")')
 580       write(*,'("(blank line required after optcomp block)")')
 581       write(*,*)
 582       stop
 583     end if
 584     if ((optcomp(1,i).lt.1).or.(optcomp(1,i).gt.3).or. &
 585         (optcomp(2,i).lt.1).or.(optcomp(2,i).gt.3).or. &
 586         (optcomp(3,i).lt.1).or.(optcomp(3,i).gt.3)) then
 587       write(*,*)
 588       write(*,'("Error(readinput): invalid optcomp : ",3I8)') optcomp
 589       write(*,*)
 590       stop
 591     end if
 592   end do
 593   write(*,*)
 594   write(*,'("Error(readinput): optical component list too long")')
 595   write(*,*)
 596   stop
 597 case('usegdft')
 598   read(50,*,err=20) usegdft
 599 case('intraband')
 600   read(50,*,err=20) intraband
 601 case('evalmin')
 602   read(50,*,err=20) evalmin
 603 case('deband')
 604   read(50,*,err=20) deband
 605   if (deband.lt.0.d0) then
 606     write(*,*)
 607     write(*,'("Error(readinput): deband < 0 : ",G18.10)') deband
 608     write(*,*)
 609     stop
 610   end if
 611 case('bfieldc')
 612   read(50,*,err=20) bfieldc
 613 case('fixspin')
 614   read(50,*,err=20) fixspin
 615 case('momfix')
 616   read(50,*,err=20) momfix
 617 case('mommtfix')
 618   do ias=1,maxspecies*maxatoms
 619     read(50,'(A80)') str
 620     if (trim(str).eq.'') goto 10
 621     read(str,*,iostat=iostat) is,ia,mommtfix(:,ia,is)
 622     if (iostat.ne.0) then
 623       write(*,*)
 624       write(*,'("Error(readinput): error reading muffin-tin fixed spin &
 625        &moments")')
 626       write(*,'("(blank line required after mommtfix block")')
 627       write(*,*)
 628       stop
 629     end if
 630   end do
 631 case('taufsm')
 632   read(50,*,err=20) taufsm
 633   if (taufsm.lt.0.d0) then
 634     write(*,*)
 635     write(*,'("Error(readinput): taufsm < 0 : ",G18.10)') taufsm
 636     write(*,*)
 637     stop
 638   end if
 639 case('autormt')
 640   read(50,*,err=20) autormt
 641 case('rmtapm')
 642   read(50,*,err=20) rmtapm
 643   if (rmtapm(1).lt.0.d0) then
 644     write(*,*)
 645     write(*,'("Error(readinput): rmtapm(1) < 0 : ",G18.10)') rmtapm(1)
 646     write(*,*)
 647     stop
 648   end if
 649   if ((rmtapm(2).le.0.d0).or.(rmtapm(2).gt.1.d0)) then
 650     write(*,*)
 651     write(*,'("Error(readinput): rmtapm(2) not in (0,1] : ",G18.10)') rmtapm(2)
 652     write(*,*)
 653     stop
 654   end if
 655 case('nosym')
 656   read(50,*,err=20) nosym
 657 case('deltaph')
 658   read(50,*,err=20) deltaph
 659 case('phwrite')
 660   read(50,*,err=20) nphwrt
 661   if (nphwrt.le.0) then
 662     write(*,*)
 663     write(*,'("Error(readinput): nphwrt <= 0 : ",I8)') nphwrt
 664     write(*,*)
 665     stop
 666   end if
 667   if (allocated(vqlwrt)) deallocate(vqlwrt)
 668   allocate(vqlwrt(3,nphwrt))
 669   do i=1,nphwrt
 670     read(50,*,err=20) vqlwrt(:,i)
 671   end do
 672 case('notes')
 673   do i=1,maxnlns
 674     read(50,'(A80)') notes(i)
 675     if (trim(notes(i)).eq.'') then
 676       notelns=i-1
 677       goto 10
 678     end if
 679   end do
 680   write(*,*)
 681   write(*,'("Error(readinput): too many note lines")')
 682   write(*,*)
 683   stop
 684 case('tforce')
 685   read(50,*,err=20) tforce
 686 case('tfibs')
 687   read(50,*,err=20) tfibs
 688 case('maxitoep')
 689   read(50,*,err=20) maxitoep
 690   if (maxitoep.lt.1) then
 691     write(*,*)
 692     write(*,'("Error(readinput): maxitoep < 1 : ",I8)') maxitoep
 693     write(*,*)
 694     stop
 695   end if
 696 case('tauoep')
 697   read(50,*,err=20) tauoep(:)
 698   if ((tauoep(1).lt.0.d0).or.(tauoep(2).lt.0.d0).or.(tauoep(3).lt.0.d0)) then
 699     write(*,*)
 700     write(*,'("Error(readinput): tauoep < 0 : ",3G18.10)') tauoep
 701     write(*,*)
 702     stop
 703   end if
 704 case('kstlist')
 705   do i=1,maxkst
 706     read(50,'(A80)') str
 707     if (trim(str).eq.'') then
 708       if (i.eq.1) then
 709         write(*,*)
 710         write(*,'("Error(readinput): empty k-point and state list")')
 711         write(*,*)
 712         stop
 713       end if
 714       nkstlist=i-1
 715       goto 10
 716     end if
 717     read(str,*,iostat=iostat) kstlist(:,i)
 718     if (iostat.ne.0) then
 719       write(*,*)
 720       write(*,'("Error(readinput): error reading k-point and state list")')
 721       write(*,'("(blank line required after kstlist block)")')
 722       write(*,*)
 723       stop
 724     end if
 725   end do
 726   write(*,*)
 727   write(*,'("Error(readinput): k-point and state list too long")')
 728   write(*,*)
 729   stop
 730 case('vklem')
 731   read(50,*,err=20) vklem
 732 case('deltaem')
 733   read(50,*,err=20) deltaem
 734 case('ndspem')
 735   read(50,*,err=20) ndspem
 736   if ((ndspem.lt.1).or.(ndspem.gt.3)) then
 737     write(*,*)
 738     write(*,'("Error(readinput): ndspem out of range : ",I8)') ndspem
 739     write(*,*)
 740     stop
 741   end if
 742 case('nosource')
 743   read(50,*,err=20) nosource
 744 case('spinsprl')
 745   read(50,*,err=20) spinsprl
 746 case('vqlss')
 747   read(50,*,err=20) vqlss
 748 case('nwrite')
 749   read(50,*,err=20) nwrite
 750 case('tevecsv')
 751   read(50,*,err=20) tevecsv
 752 case('lda+u')
 753   read(50,*) ldapu
 754   do is=1,maxspecies
 755     read(50,'(A80)') str
 756     if (trim(str).eq.'') goto 10
 757     read(str,*,iostat=iostat) js,l,t1,t2
 758     if (iostat.ne.0) then
 759       write(*,*)
 760       write(*,'("Error(readinput): error reading LDA+U parameters")')
 761       write(*,'("(blank line required after lda+u block)")')
 762       write(*,*)
 763       stop
 764     end if
 765     if ((js.le.0).or.(js.ge.maxspecies)) then
 766       write(*,*)
 767       write(*,'("Error(readinput): invalid species number in lda+u block : ",&
 768        &I8)') js
 769       write(*,*)
 770       stop
 771     end if
 772     if (l.gt.lmaxlu) then
 773       write(*,*)
 774       write(*,'("Error(readinput): l > lmaxlu in lda+u block : ",2I8)') l,lmaxlu
 775       write(*,*)
 776       stop
 777     end if
 778     llu(js)=l
 779     ujlu(1,js)=t1
 780     ujlu(2,js)=t2
 781   end do
 782 case('rdmxctype')
 783   read(50,*,err=20) rdmxctype
 784 case('rdmmaxscl')
 785   read(50,*,err=20) rdmmaxscl
 786   if (rdmmaxscl.lt.0) then
 787     write(*,*)
 788     write(*,'("Error(readinput): rdmmaxscl < 0 : ",I8)') rdmmaxscl
 789     write(*,*)
 790   end if
 791 case('maxitn')
 792   read(50,*,err=20) maxitn
 793 case('maxitc')
 794   read(50,*,err=20) maxitc
 795 case('taurdmn')
 796   read(50,*,err=20) taurdmn
 797   if (taurdmn.lt.0.d0) then
 798     write(*,*)
 799     write(*,'("Error(readinput): taurdmn < 0 : ",G18.10)') taurdmn
 800     write(*,*)
 801     stop
 802   end if
 803 case('taurdmc')
 804   read(50,*,err=20) taurdmc
 805   if (taurdmc.lt.0.d0) then
 806     write(*,*)
 807     write(*,'("Error(readinput): taurdmc < 0 : ",G18.10)') taurdmc
 808     write(*,*)
 809     stop
 810   end if
 811 case('rdmalpha')
 812   read(50,*,err=20) rdmalpha
 813   if ((rdmalpha.le.0.d0).or.(rdmalpha.ge.1.d0)) then
 814     write(*,*)
 815     write(*,'("Error(readinput): rdmalpha not in (0,1) : ",G18.10)') rdmalpha
 816     write(*,*)
 817     stop
 818   end if
 819 case('reducebf')
 820   read(50,*,err=20) reducebf
 821   if ((reducebf.lt.0.d0).or.(reducebf.gt.1.d0)) then
 822     write(*,*)
 823     write(*,'("Error(readinput): reducebf not in [0,1] : ",G18.10)') reducebf
 824     write(*,*)
 825     stop
 826   end if
 827 case('')
 828   goto 10
 829 case default
 830   write(*,*)
 831   write(*,'("Error(readinput): invalid block name : ",A)') trim(bname)
 832   write(*,*)
 833   stop
 834 end select
 835 goto 10
 836 20 continue
 837 write(*,*)
 838 write(*,'("Error(readinput): error reading from exciting.in")')
 839 write(*,'("Problem occurred in ''",A,"'' block")') trim(bname)
 840 write(*,*)
 841 stop
 842 30 continue
 843 close(50)
 844 ! scale the lattice vectors (scaling not referenced again in code)
 845 avec(:,1)=sc1*avec(:,1)
 846 avec(:,2)=sc2*avec(:,2)
 847 avec(:,3)=sc3*avec(:,3)
 848 avec(:,:)=sc*avec(:,:)
 849 ! check if system is an isolated molecule
 850 if (molecule) then
 851 ! set up cubic unit cell with vacuum region around molecule
 852   avec(:,:)=0.d0
 853   do is=1,nspecies
 854     do ia=1,natoms(is)
 855       do js=1,nspecies
 856         do ja=1,natoms(is)
 857           do i=1,3
 858             t1=abs(atposl(i,ia,is)-atposl(i,ja,js))
 859             if (t1.gt.avec(i,i)) avec(i,i)=t1
 860           end do
 861         end do
 862       end do
 863     end do
 864   end do
 865   do i=1,3
 866     avec(i,i)=avec(i,i)+vacuum
 867   end do
 868 ! convert atomic positions from Cartesian to lattice coordinates
 869   call r3minv(avec,ainv)
 870   do is=1,nspecies
 871     do ia=1,natoms(is)
 872       call r3mv(ainv,atposl(1,ia,is),v)
 873       atposl(:,ia,is)=v(:)
 874     end do
 875   end do
 876 end if
 877
 878 !---------------------------------------------!
 879 !     read from atomic species data files     !
 880 !---------------------------------------------!
 881 do is=1,nspecies
 882   open(50,file=trim(sppath)//trim(spfname(is)),action='READ',status='OLD', &
 883    form='FORMATTED',iostat=iostat)
 884   if (iostat.ne.0) then
 885     write(*,*)
 886     write(*,'("Error(readinput): error opening species file ",A)') &
 887      trim(sppath)//trim(spfname(is))
 888     write(*,*)
 889     stop
 890   end if
 891   read(50,*) spsymb(is)
 892   read(50,*) spname(is)
 893   read(50,*) spzn(is)
 894   read(50,*) spmass(is)
 895   read(50,*) sprmin(is),rmt(is),sprmax(is),nrmt(is)
 896   if (sprmin(is).le.0.d0) then
 897     write(*,*)
 898     write(*,'("Error(readinput): sprmin <= 0 : ",G18.10)') sprmin(is)
 899     write(*,'(" for species ",I4)') is
 900     write(*,*)
 901     stop
 902   end if
 903   if (rmt(is).le.sprmin(is)) then
 904     write(*,*)
 905     write(*,'("Error(readinput): rmt <= sprmin : ",2G18.10)') rmt(is),sprmin(is)
 906     write(*,'(" for species ",I4)') is
 907     write(*,*)
 908     stop
 909   end if
 910   if (sprmax(is).lt.rmt(is)) then
 911     write(*,*)
 912     write(*,'("Error(readinput): sprmax < rmt : ",2G18.10)') sprmax(is),rmt(is)
 913     write(*,*)
 914     stop
 915   end if
 916   if (nrmt(is).lt.20) then
 917     write(*,*)
 918     write(*,'("Error(readinput): nrmt too small : ",I8)') nrmt(is)
 919     write(*,'(" for species ",I4)') is
 920     write(*,*)
 921     stop
 922   end if
 923   read(50,*) spnst(is)
 924   if (spnst(is).le.0) then
 925     write(*,*)
 926     write(*,'("Error(readinput): invalid spnst : ",I8)') spnst(is)
 927     write(*,'(" for species ",I4)') is
 928     write(*,*)
 929     stop
 930   end if
 931   if (spnst(is).gt.maxspst) then
 932     write(*,*)
 933     write(*,'("Error(readinput): too many states for species ",I8)') is
 934     write(*,*)
 935     stop
 936   end if
 937   do ist=1,spnst(is)
 938     read(50,*) spn(ist,is),spl(ist,is),spk(ist,is),spocc(ist,is),spcore(ist,is)
 939     if (spn(ist,is).lt.1) then
 940       write(*,*)
 941       write(*,'("Error(readinput): spn < 1 : ",I8)') spn(ist,is)
 942       write(*,'(" for species ",I4)') is
 943       write(*,'(" and state ",I4)') ist
 944       write(*,*)
 945       stop
 946     end if
 947     if (spl(ist,is).lt.0) then
 948       write(*,*)
 949       write(*,'("Error(readinput): spl < 0 : ",I8)') spl(ist,is)
 950       write(*,'(" for species ",I4)') is
 951       write(*,'(" and state ",I4)') ist
 952       write(*,*)
 953       stop
 954     end if
 955     if (spk(ist,is).lt.1) then
 956       write(*,*)
 957       write(*,'("Error(readinput): spk < 1 : ",I8)') spk(ist,is)
 958       write(*,'(" for species ",I4)') is
 959       write(*,'(" and state ",I4)') ist
 960       write(*,*)
 961       stop
 962     end if
 963     if (spocc(ist,is).lt.0.d0) then
 964       write(*,*)
 965       write(*,'("Error(readinput): spocc < 0 : ",G18.10)') spocc(ist,is)
 966       write(*,'(" for species ",I4)') is
 967       write(*,'(" and state ",I4)') ist
 968       write(*,*)
 969       stop
 970     end if
 971   end do
 972   read(50,*) apword(0,is)
 973   if (apword(0,is).le.0) then
 974     write(*,*)
 975     write(*,'("Error(readinput): apword <= 0 : ",I8)') apword(0,is)
 976     write(*,'(" for species ",I4)') is
 977     write(*,*)
 978     stop
 979   end if
 980   if (apword(0,is).gt.maxapword) then
 981     write(*,*)
 982     write(*,'("Error(readinput): apword too large : ",I8)') apword(0,is)
 983     write(*,'(" for species ",I4)') is
 984     write(*,'("Adjust maxapword in modmain and recompile code")')
 985     write(*,*)
 986     stop
 987   end if
 988 ! set the APW orders for l>0
 989   apword(1:lmaxapw,is)=apword(0,is)
 990   do io=1,apword(0,is)
 991     read(50,*) apwe0(io,0,is),apwdm(io,0,is),apwve(io,0,is)
 992     if (apwdm(io,0,is).lt.0) then
 993       write(*,*)
 994       write(*,'("Error(readinput): apwdm < 0 : ",I8)') apwdm(io,0,is)
 995       write(*,'(" for species ",I4)') is
 996       write(*,'(" and order ",I4)') io
 997       write(*,*)
 998       stop
 999     end if
1000 ! set the APW linearisation energies, derivative orders and variability for l>0
1001     apwe0(io,1:lmaxapw,is)=apwe0(io,0,is)
1002     apwdm(io,1:lmaxapw,is)=apwdm(io,0,is)
1003     apwve(io,1:lmaxapw,is)=apwve(io,0,is)
1004   end do
1005   read(50,*) nlx
1006   if (nlx.lt.0) then
1007     write(*,*)
1008     write(*,'("Error(readinput): nlx < 0 : ",I8)') nlx
1009     write(*,'(" for species ",I4)') is
1010     write(*,*)
1011     stop
1012   end if
1013   do ilx=1,nlx
1014     read(50,*) lx,io
1015     if (lx.lt.0) then
1016       write(*,*)
1017       write(*,'("Error(readinput): lx < 0 : ",I8)') lx
1018       write(*,'(" for species ",I4)') is
1019       write(*,'(" and exception number ",I4)') ilx
1020       write(*,*)
1021       stop
1022     end if
1023     if (lx.gt.lmaxapw) then
1024       write(*,*)
1025       write(*,'("Error(readinput): lx > lmaxapw : ",I8)') lx
1026       write(*,'(" for species ",I4)') is
1027       write(*,'(" and exception number ",I4)') ilx
1028       write(*,*)
1029       stop
1030     end if
1031     apword(lx,is)=io
1032     if (apword(lx,is).le.0) then
1033       write(*,*)
1034       write(*,'("Error(readinput): apword <= 0 : ",I8)') apword(lx,is)
1035       write(*,'(" for species ",I4)') is
1036       write(*,'(" and exception number ",I4)') ilx
1037       write(*,*)
1038       stop
1039     end if
1040     if (apword(lx,is).gt.maxapword) then
1041       write(*,*)
1042       write(*,'("Error(readinput): apword too large : ",I8)') apword(lx,is)
1043       write(*,'(" for species ",I4)') is
1044       write(*,'(" and exception number ",I4)') ilx
1045       write(*,'("Adjust maxapword in modmain and recompile code")')
1046       write(*,*)
1047       stop
1048     end if
1049     do io=1,apword(lx,is)
1050       read(50,*) apwe0(io,lx,is),apwdm(io,lx,is),apwve(io,lx,is)
1051       if (apwdm(io,lx,is).lt.0) then
1052         write(*,*)
1053         write(*,'("Error(readinput): apwdm < 0 : ",I8)') apwdm(io,lx,is)
1054         write(*,'(" for species ",I4)') is
1055         write(*,'(" exception number ",I4)') ilx
1056         write(*,'(" and order ",I4)') io
1057         write(*,*)
1058         stop
1059       end if
1060     end do
1061   end do
1062   read(50,*) nlorb(is)
1063   if (nlorb(is).lt.0) then
1064     write(*,*)
1065     write(*,'("Error(readinput): nlorb < 0 : ",I8)') nlorb(is)
1066     write(*,'(" for species ",I4)') is
1067     write(*,*)
1068     stop
1069   end if
1070   if (nlorb(is).gt.maxlorb) then
1071     write(*,*)
1072     write(*,'("Error(readinput): nlorb too large : ",I8)') nlorb(is)
1073     write(*,'(" for species ",I4)') is
1074     write(*,'("Adjust maxlorb in modmain and recompile code")')
1075     write(*,*)
1076     stop
1077   end if
1078   do ilo=1,nlorb(is)
1079     read(50,*) lorbl(ilo,is),lorbord(ilo,is)
1080     if (lorbl(ilo,is).lt.0) then
1081       write(*,*)
1082       write(*,'("Error(readinput): lorbl < 0 : ",I8)') lorbl(ilo,is)
1083       write(*,'(" for species ",I4)') is
1084       write(*,'(" and local-orbital ",I4)') ilo
1085       write(*,*)
1086       stop
1087     end if
1088     if (lorbl(ilo,is).gt.lmaxmat) then
1089       write(*,*)
1090       write(*,'("Error(readinput): lorbl > lmaxmat : ",2I8)') lorbl(ilo,is), &
1091        lmaxmat
1092       write(*,'(" for species ",I4)') is
1093       write(*,'(" and local-orbital ",I4)') ilo
1094       write(*,*)
1095       stop
1096     end if
1097     if (lorbord(ilo,is).lt.2) then
1098       write(*,*)
1099       write(*,'("Error(readinput): lorbord < 2 : ",I8)') lorbord(ilo,is)
1100       write(*,'(" for species ",I4)') is
1101       write(*,'(" and local-orbital ",I4)') ilo
1102       write(*,*)
1103       stop
1104     end if
1105     if (lorbord(ilo,is).gt.maxlorbord) then
1106       write(*,*)
1107       write(*,'("Error(readinput): lorbord too large : ",I8)') lorbord(ilo,is)
1108       write(*,'(" for species ",I4)') is
1109       write(*,'(" and local-orbital ",I4)') ilo
1110       write(*,'("Adjust maxlorbord in modmain and recompile code")')
1111       write(*,*)
1112       stop
1113     end if
1114     do io=1,lorbord(ilo,is)
1115       read(50,*) lorbe0(io,ilo,is),lorbdm(io,ilo,is),lorbve(io,ilo,is)
1116       if (lorbdm(io,ilo,is).lt.0) then
1117         write(*,*)
1118         write(*,'("Error(readinput): lorbdm < 0 : ",I8)') lorbdm(io,ilo,is)
1119         write(*,'(" for species ",I4)') is
1120         write(*,'(" local-orbital ",I4)') ilo
1121         write(*,'(" and order ",I4)') io
1122         write(*,*)
1123         stop
1124       end if
1125     end do
1126   end do
1127   close(50)
1128 end do
1129 end subroutine
1130 !EOC