1 2! Copyright (C) 2002-2009 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 6module modmain 7 8!----------------------------! 9! lattice parameters ! 10!----------------------------! 11! lattice vectors stored column-wise 12real(8) avec(3,3) 13! magnitude of random displacements added to lattice vectors 14real(8) rndavec 15! inverse of lattice vector matrix 16real(8) ainv(3,3) 17! reciprocal lattice vectors 18real(8) bvec(3,3) 19! inverse of reciprocal lattice vector matrix 20real(8) binv(3,3) 21! unit cell volume 22real(8) omega 23! Brillouin zone volume 24real(8) omegabz 25! any vector with length less than epslat is considered zero 26real(8) epslat 27 28!--------------------------! 29! atomic variables ! 30!--------------------------! 31! maximum allowed species 32integer, parameter :: maxspecies=8 33! maximum allowed atoms per species 34integer, parameter :: maxatoms=200 35! number of species 36integer nspecies 37! number of atoms for each species 38integer natoms(maxspecies) 39! maximum number of atoms over all the species 40integer natmmax 41! total number of atoms 42integer natmtot 43! index to atoms and species 44integer idxas(maxatoms,maxspecies) 45! inverse atoms and species indices 46integer idxis(maxatoms*maxspecies) 47integer idxia(maxatoms*maxspecies) 48! molecule is .true. is the system is an isolated molecule 49logical molecule 50! primcell is .true. if primitive unit cell is to be found automatically 51logical primcell 52! atomic positions in lattice coordinates 53real(8) atposl(3,maxatoms,maxspecies) 54! atomic positions in Cartesian coordinates 55real(8) atposc(3,maxatoms,maxspecies) 56! magnitude of random displacements added to the atomic positions 57real(8) rndatposc 58! tdatpos is .true. if small amplitude displacements are to be included when 59! calculating the Coulomb potential 60logical :: tdatpos=.false. 61! trddatpos is .true. if the small amplitude displacements and velocities are to 62! be read from file 63logical :: trddatpos=.false. 64! small amplitude atomic displacements and velocities 65real(8) datposc(3,0:1,maxatoms,maxspecies) 66 67!----------------------------------! 68! atomic species variables ! 69!----------------------------------! 70! species files path 71character(256) sppath 72! species filenames 73character(256) spfname(maxspecies) 74! species name 75character(256) spname(maxspecies) 76! species symbol 77character(256) spsymb(maxspecies) 78! species nuclear charge 79real(8) spzn(maxspecies) 80! ptnucl is .true. if the nuclei are to be treated as point charges, if .false. 81! the nuclei have a finite spherical distribution 82logical ptnucl 83! nuclear radius 84real(8) rnucl(maxspecies) 85! nuclear volume 86real(8) volnucl(maxspecies) 87! number of radial mesh points to nuclear radius 88integer nrnucl(maxspecies) 89! nuclear Coulomb potential 90real(8), allocatable :: vcln(:,:) 91! species electronic charge 92real(8) spze(maxspecies) 93! species mass 94real(8) spmass(maxspecies) 95! smallest radial point for each species 96real(8) rminsp(maxspecies) 97! effective infinity for species 98real(8) rmaxsp(maxspecies) 99! number of radial points to effective infinity for each species 100integer nrsp(maxspecies) 101! maximum nrsp over all the species 102integer nrspmax 103! maximum allowed states for each species 104integer, parameter :: maxstsp=40 105! number of states for each species 106integer nstsp(maxspecies) 107! maximum nstsp over all the species 108integer nstspmax 109! core-valence cut-off energy for species file generation 110real(8) ecvcut 111! semi-core-valence cut-off energy for species file generation 112real(8) esccut 113! state principle quantum number for each species 114integer nsp(maxstsp,maxspecies) 115! state l value for each species 116integer lsp(maxstsp,maxspecies) 117! state k value for each species 118integer ksp(maxstsp,maxspecies) 119! spcore is .true. if species state is core 120logical spcore(maxstsp,maxspecies) 121! total number of core states 122integer nstcr 123! state eigenvalue for each species 124real(8) evalsp(maxstsp,maxspecies) 125! state occupancy for each species 126real(8) occsp(maxstsp,maxspecies) 127! species radial mesh to effective infinity 128real(8), allocatable :: rsp(:,:) 129! species charge density 130real(8), allocatable :: rhosp(:,:) 131! species self-consistent potential 132real(8), allocatable :: vrsp(:,:) 133! exchange-correlation type for atomic species (the converged ground-state of 134! the crystal does not depend on this choice) 135integer xctsp(3) 136 137!---------------------------------------------------------------! 138! muffin-tin radial mesh and angular momentum variables ! 139!---------------------------------------------------------------! 140! scale factor for number of muffin-tin points 141real(8) nrmtscf 142! number of muffin-tin radial points for each species 143integer nrmt(maxspecies) 144! maximum nrmt over all the species 145integer nrmtmax 146! optional default muffin-tin radius for all atoms 147real(8) rmtall 148! minimum allowed distance between muffin-tin surfaces 149real(8) rmtdelta 150! muffin-tin radii 151real(8) rmt(maxspecies) 152! total muffin-tin volume 153real(8) omegamt 154! radial step length for coarse mesh 155integer lradstp 156! number of coarse radial mesh points 157integer nrcmt(maxspecies) 158! maximum nrcmt over all the species 159integer nrcmtmax 160! coarse muffin-tin radial mesh 161real(8), allocatable :: rcmt(:,:) 162! r^l on fine radial mesh 163real(8), allocatable :: rlmt(:,:,:) 164! r^l on coarse radial mesh 165real(8), allocatable :: rlcmt(:,:,:) 166! weights for spline integration on fine radial mesh 167real(8), allocatable :: wrmt(:,:) 168! weights for spline partial integration on fine radial mesh 169real(8), allocatable :: wprmt(:,:,:) 170! weights for spline coefficients on fine radial mesh 171real(8), allocatable :: wcrmt(:,:,:) 172! weights for spline integration on coarse radial mesh 173real(8), allocatable :: wrcmt(:,:) 174! weights for spline partial integration on coarse radial mesh 175real(8), allocatable :: wprcmt(:,:,:) 176! weights for spline coefficients on coarse radial mesh 177real(8), allocatable :: wcrcmt(:,:,:) 178! maximum allowable angular momentum for augmented plane waves 179integer, parameter :: maxlapw=30 180! maximum angular momentum for augmented plane waves 181integer lmaxapw 182! (lmaxapw+1)^2 183integer lmmaxapw 184! maximum angular momentum on the outer part of the muffin-tin 185integer lmaxo 186! (lmaxo+1)^2 187integer lmmaxo 188! maximum angular momentum on the inner part of the muffin-tin 189integer lmaxi 190! (lmaxi+1)^2 191integer lmmaxi 192! fraction of muffin-tin radius which constitutes the inner part 193real(8) fracinr 194! number of fine/coarse radial points on the inner part of the muffin-tin 195integer nrmti(maxspecies),nrcmti(maxspecies) 196! number of fine/coarse points in packed muffin-tins 197integer npmti(maxspecies),npmt(maxspecies) 198integer npcmti(maxspecies),npcmt(maxspecies) 199! maximum number of points over all packed muffin-tins 200integer npmtmax,npcmtmax 201 202!--------------------------------! 203! spin related variables ! 204!--------------------------------! 205! spinpol is .true. for spin-polarised calculations 206logical spinpol 207! spinorb is .true. for spin-orbit coupling 208logical spinorb 209! scale factor of spin-orbit coupling term in Hamiltonian 210real(8) socscf 211! dimension of magnetisation and magnetic vector fields (1 or 3) 212integer ndmag 213! ncmag is .true. if the magnetisation is non-collinear, i.e. when ndmag = 3 214logical ncmag 215! if cmagz is .true. then collinear magnetism along the z-axis is enforced 216logical cmagz 217! spcpl is .true. if the up and down spins are coupled 218logical spcpl 219! fixed spin moment type 220! 0 : none 221! 1 (-1) : total moment (direction) 222! 2 (-2) : individual muffin-tin moments (direction) 223! 3 (-3) : total and muffin-tin moments (direction) 224integer fsmtype 225! fixed total spin magnetic moment 226real(8) momfix(3) 227! fixed spin moment global effective field in Cartesian coordinates 228real(8) bfsmc(3) 229! muffin-tin fixed spin moments 230real(8) mommtfix(3,maxatoms,maxspecies) 231! muffin-tin fixed spin moment effective fields in Cartesian coordinates 232real(8), allocatable :: bfsmcmt(:,:) 233! fixed spin moment field step size 234real(8) taufsm 235! second-variational spinor dimension (1 or 2) 236integer nspinor 237! global external magnetic field in Cartesian coordinates 238real(8) bfieldc(3) 239! initial field 240real(8) bfieldc0(3) 241! external magnetic field in each muffin-tin in Cartesian coordinates 242real(8) bfcmt(3,maxatoms,maxspecies) 243! initial field as read in from input file 244real(8) bfcmt0(3,maxatoms,maxspecies) 245! magnitude of random vectors added to muffin-tin fields 246real(8) rndbfcmt 247! external magnetic fields are multiplied by reducebf after each s.c. loop 248real(8) reducebf 249! spinsprl is .true. if a spin-spiral is to be calculated 250logical spinsprl 251! ssdph is .true. if the muffin-tin spin-spiral magnetisation is de-phased 252logical ssdph 253! spin-spiral phase factor for each atom 254complex(8), allocatable :: zqss(:) 255! number of spin-dependent first-variational functions per state 256integer nspnfv 257! map from second- to first-variational spin index 258integer jspnfv(2) 259! spin-spiral q-vector in lattice coordinates 260real(8) vqlss(3) 261! spin-spiral q-vector in Cartesian coordinates 262real(8) vqcss(3) 263! current q-point in spin-spiral supercell calculation 264integer iqss 265! number of primitive unit cells in spin-spiral supercell 266integer nscss 267! number of fixed spin direction points on the sphere for finding the magnetic 268! anisotropy energy (MAE) 269integer npmae0,npmae 270! (theta,phi) coordinates for each MAE direction 271real(8), allocatable :: tpmae(:,:) 272 273!---------------------------------------------! 274! electric field and vector potential ! 275!---------------------------------------------! 276! tefield is .true. if a polarising constant electric field is applied 277logical tefield 278! electric field vector in Cartesian coordinates 279real(8) efieldc(3) 280! electric field vector in lattice coordinates 281real(8) efieldl(3) 282! tafield is .true. if a constant vector potential is applied 283logical tafield 284! vector potential A-field which couples to paramagnetic current 285real(8) afieldc(3) 286! A-field in lattice coordinates 287real(8) afieldl(3) 288 289!----------------------------! 290! symmetry variables ! 291!----------------------------! 292! type of symmetry allowed for the crystal 293! 0 : only the identity element is used 294! 1 : full symmetry group is used 295! 2 : only symmorphic symmetries are allowed 296integer symtype 297! number of Bravais lattice point group symmetries 298integer nsymlat 299! Bravais lattice point group symmetries 300integer symlat(3,3,48) 301! determinants of lattice symmetry matrices (1 or -1) 302integer symlatd(48) 303! index to inverses of the lattice symmetries 304integer isymlat(48) 305! lattice point group symmetries in Cartesian coordinates 306real(8) symlatc(3,3,48) 307! tshift is .true. if atomic basis is allowed to be shifted 308logical tshift 309! tsyminv is .true. if the crystal has inversion symmetry 310logical tsyminv 311! maximum of symmetries allowed 312integer, parameter :: maxsymcrys=192 313! number of crystal symmetries 314integer nsymcrys 315! crystal symmetry translation vector in lattice and Cartesian coordinates 316real(8) vtlsymc(3,maxsymcrys) 317real(8) vtcsymc(3,maxsymcrys) 318! tv0symc is .true. if the translation vector is zero 319logical tv0symc(maxsymcrys) 320! spatial rotation element in lattice point group for each crystal symmetry 321integer lsplsymc(maxsymcrys) 322! global spin rotation element in lattice point group for each crystal symmetry 323integer lspnsymc(maxsymcrys) 324! equivalent atom index for each crystal symmetry 325integer, allocatable :: ieqatom(:,:,:) 326! eqatoms(ia,ja,is) is .true. if atoms ia and ja are equivalent 327logical, allocatable :: eqatoms(:,:,:) 328! number of site symmetries 329integer, allocatable :: nsymsite(:) 330! site symmetry spatial rotation element in lattice point group 331integer, allocatable :: lsplsyms(:,:) 332! site symmetry global spin rotation element in lattice point group 333integer, allocatable :: lspnsyms(:,:) 334 335!----------------------------! 336! G-vector variables ! 337!----------------------------! 338! G-vector cut-off for interstitial potential and density 339real(8) gmaxvr 340! G-vector grid sizes 341integer ngridg(3) 342! G-vector grid sizes for coarse grid (G < 2*gkmax) 343integer ngdgc(3) 344! total number of G-vectors 345integer ngtot 346! total number of G-vectors for coarse grid (G < 2*gkmax) 347integer ngtc 348! integer grid intervals for each direction 349integer intgv(2,3) 350! number of G-vectors with G < gmaxvr 351integer ngvec 352! number of G-vectors for coarse grid (G < 2*gkmax) 353integer ngvc 354! G-vector integer coordinates (i1,i2,i3) 355integer, allocatable :: ivg(:,:) 356! map from (i1,i2,i3) to G-vector index 357integer, allocatable :: ivgig(:,:,:) 358! map from G-vector index to FFT array 359integer, allocatable :: igfft(:) 360! map from G-vector index to FFT array for coarse grid (G < 2*gkmax) 361integer, allocatable :: igfc(:) 362! G-vectors in Cartesian coordinates 363real(8), allocatable :: vgc(:,:) 364! length of G-vectors 365real(8), allocatable :: gc(:) 366! Coulomb Green's function in G-space = 4 pi / G^2 367real(8), allocatable :: gclg(:) 368! spherical Bessel functions j_l(|G|R_mt) 369real(8), allocatable :: jlgrmt(:,:,:) 370! spherical harmonics of the G-vectors 371complex(8), allocatable :: ylmg(:,:) 372! structure factors for the G-vectors 373complex(8), allocatable :: sfacg(:,:) 374! smooth step function form factors for all species and G-vectors 375real(8), allocatable :: ffacg(:,:) 376! characteristic function in G-space: 0 inside the muffin-tins and 1 outside 377complex(8), allocatable :: cfunig(:) 378! characteristic function in real-space: 0 inside the muffin-tins and 1 outside 379real(8), allocatable :: cfunir(:) 380! characteristic function in real-space for coarse grid (G < 2*gkmax) 381real(8), allocatable :: cfrc(:) 382 383!---------------------------! 384! k-point variables ! 385!---------------------------! 386! autokpt is .true. if the k-point set is determined automatically 387logical autokpt 388! radius of sphere used to determine k-point density when autokpt is .true. 389real(8) radkpt 390! k-point grid sizes 391integer ngridk(3) 392! k-point offset 393real(8) vkloff(3) 394! type of reduction to perform on k-point set 395! 0 : no reduction 396! 1 : reduce with full crystal symmetry group 397! 2 : reduce with symmorphic symmetries only 398integer reducek 399! number of point group symmetries used for k-point reduction 400integer nsymkpt 401! point group symmetry matrices used for k-point reduction 402integer symkpt(3,3,48) 403! total number of reduced k-points 404integer nkpt 405! total number of non-reduced k-points 406integer nkptnr 407! locations of k-points on integer grid 408integer, allocatable :: ivk(:,:) 409! map from integer grid to reduced k-point index 410integer, allocatable :: ivkik(:,:,:) 411! map from integer grid to non-reduced k-point index 412integer, allocatable :: ivkiknr(:,:,:) 413! k-points in lattice coordinates 414real(8), allocatable :: vkl(:,:) 415! k-points in Cartesian coordinates 416real(8), allocatable :: vkc(:,:) 417! reduced k-point weights 418real(8), allocatable :: wkpt(:) 419! weight of each non-reduced k-point 420real(8) wkptnr 421! k-point at which to determine effective mass tensor 422real(8) vklem(3) 423! displacement size for computing the effective mass tensor 424real(8) deltaem 425! number of displacements in each direction 426integer ndspem 427! number of k-points subdivision used for calculating the polarisation phase 428integer nkspolar 429 430!------------------------------! 431! G+k-vector variables ! 432!------------------------------! 433! species for which the muffin-tin radius will be used for calculating gkmax 434integer isgkmax 435! smallest muffin-tin radius times gkmax 436real(8) rgkmax 437! maximum |G+k| cut-off for APW functions 438real(8) gkmax 439! number of G+k-vectors for augmented plane waves 440integer, allocatable :: ngk(:,:) 441! maximum number of G+k-vectors over all k-points 442integer ngkmax 443! index from G+k-vectors to G-vectors 444integer, allocatable :: igkig(:,:,:) 445! G+k-vectors in lattice coordinates 446real(8), allocatable :: vgkl(:,:,:,:) 447! G+k-vectors in Cartesian coordinates 448real(8), allocatable :: vgkc(:,:,:,:) 449! length of G+k-vectors 450real(8), allocatable :: gkc(:,:,:) 451! structure factors for the G+k-vectors 452complex(8), allocatable :: sfacgk(:,:,:,:) 453 454!---------------------------! 455! q-point variables ! 456!---------------------------! 457! q-point grid sizes 458integer ngridq(3) 459! integer grid intervals for the q-points 460integer intq(2,3) 461! type of reduction to perform on q-point set (see reducek) 462integer reduceq 463! number of point group symmetries used for q-point reduction 464integer nsymqpt 465! point group symmetry matrices used for q-point reduction 466integer symqpt(3,3,48) 467! total number of reduced q-points 468integer nqpt 469! total number of non-reduced q-points 470integer nqptnr 471! locations of q-points on integer grid 472integer, allocatable :: ivq(:,:) 473! map from integer grid to reduced index 474integer, allocatable :: ivqiq(:,:,:) 475! map from integer grid to non-reduced index 476integer, allocatable :: ivqiqnr(:,:,:) 477! map from q-vector index to complex-complex FFT array 478integer, allocatable :: iqfft(:) 479! number of complex FFT elements for real-complex transforms 480integer nfqrz 481! map from q-point index to real-complex FFT index 482integer, allocatable :: ifqrz(:) 483! map from real-complex FFT index to q-point index 484integer, allocatable :: iqrzf(:) 485! q-points in lattice coordinates 486real(8), allocatable :: vql(:,:) 487! q-points in Cartesian coordinates 488real(8), allocatable :: vqc(:,:) 489! q-point weights 490real(8), allocatable :: wqpt(:) 491! weight for each non-reduced q-point 492real(8) wqptnr 493! regularised Coulomb Green's function in q-space 494real(8), allocatable :: gclq(:) 495! if t0gclq0 is .true. then the Coulomb Green's function at q = 0 is set to zero 496logical t0gclq0 497 498!-----------------------------------------------------! 499! spherical harmonic transform (SHT) matrices ! 500!-----------------------------------------------------! 501! trotsht is .true. if the spherical cover used for the SHT is to be rotated 502logical :: trotsht=.false. 503! spherical cover rotation matrix 504real(8) rotsht(3,3) 505! real backward SHT matrix for lmaxi 506real(8), allocatable :: rbshti(:,:) 507! real forward SHT matrix for lmaxi 508real(8), allocatable :: rfshti(:,:) 509! real backward SHT matrix for lmaxo 510real(8), allocatable :: rbshto(:,:) 511! real forward SHT matrix for lmaxo 512real(8), allocatable :: rfshto(:,:) 513! complex backward SHT matrix for lmaxi 514complex(8), allocatable :: zbshti(:,:) 515! complex forward SHT matrix for lmaxi 516complex(8), allocatable :: zfshti(:,:) 517! complex backward SHT matrix for lmaxo 518complex(8), allocatable :: zbshto(:,:) 519! complex forward SHT matrix for lmaxo 520complex(8), allocatable :: zfshto(:,:) 521 522!---------------------------------------------------------------! 523! density, potential and exchange-correlation variables ! 524!---------------------------------------------------------------! 525! exchange-correlation functional type 526integer xctype(3) 527! exchange-correlation functional description 528character(512) xcdescr 529! exchange-correlation functional spin requirement 530integer xcspin 531! exchange-correlation functional density gradient requirement 532integer xcgrad 533! small constant used to stabilise non-collinear GGA 534real(8) dncgga 535! muffin-tin and interstitial charge density 536real(8), allocatable :: rhomt(:,:),rhoir(:) 537! trhonorm is .true. if the density is to be normalised after every iteration 538logical trhonorm 539! muffin-tin and interstitial magnetisation vector field 540real(8), allocatable :: magmt(:,:,:),magir(:,:) 541! tjr is .true. if the current density j(r) is to be calculated 542logical tjr 543! muffin-tin and interstitial gauge-invariant current density vector field 544real(8), allocatable :: jrmt(:,:,:),jrir(:,:) 545! amount of smoothing to be applied to the exchange-correlation potentials and 546! magnetic field 547integer msmooth 548! muffin-tin and interstitial Coulomb potential 549real(8), allocatable :: vclmt(:,:),vclir(:) 550! Poisson solver pseudocharge density constant 551integer npsd 552! lmaxo+npsd+1 553integer lnpsd 554! muffin-tin and interstitial exchange energy density 555real(8), allocatable :: exmt(:,:),exir(:) 556! muffin-tin and interstitial correlation energy density 557real(8), allocatable :: ecmt(:,:),ecir(:) 558! muffin-tin and interstitial exchange-correlation potential 559real(8), allocatable :: vxcmt(:,:),vxcir(:) 560! muffin-tin and interstitial exchange-correlation magnetic field 561real(8), allocatable :: bxcmt(:,:,:),bxcir(:,:) 562! muffin-tin and interstitial magnetic dipole field 563real(8), allocatable :: bdmt(:,:,:),bdir(:,:) 564! tbdip is .true. if the spin and current dipole fields are to be added to the 565! Kohn-Sham magnetic field 566logical tbdip 567! combined target array for vsmt, vsir, bsmt and bsir 568real(8), allocatable, target :: vsbs(:) 569! muffin-tin and interstitial Kohn-Sham effective potential 570real(8), pointer :: vsmt(:,:),vsir(:) 571! muffin-tin Kohn-Sham effective magnetic field in spherical coordinates and on 572! a coarse radial mesh 573real(8), pointer :: bsmt(:,:,:) 574! interstitial Kohn-Sham effective magnetic field 575real(8), pointer :: bsir(:,:) 576! G-space interstitial Kohn-Sham effective potential 577complex(8), allocatable :: vsig(:) 578! nosource is .true. if the field is to be made source-free 579logical nosource 580! tssxc is .true. if scaled spin exchange-correlation is to be used 581logical tssxc 582! spin exchange-correlation scaling factor 583real(8) sxcscf 584! spin-orbit coupling radial function 585real(8), allocatable :: socfr(:,:) 586! kinetic energy density 587real(8), allocatable :: taumt(:,:,:),tauir(:,:) 588! core kinetic energy density 589real(8), allocatable :: taucr(:,:,:) 590! taudft is .true. if meta-GGA is to be treated as a meta-GGA functional 591logical taudft 592! meta-GGA exchange-correlation potential 593real(8), allocatable :: wxcmt(:,:),wxcir(:) 594! meta-GGA Kohn-Sham potential 595real(8), allocatable :: wsmt(:,:),wsir(:) 596! Tran-Blaha '09 constant c [Phys. Rev. Lett. 102, 226401 (2009)] 597real(8) c_tb09 598! tc_tb09 is .true. if the Tran-Blaha constant has been read in 599logical tc_tb09 600! if trdstate is .true. the density and potential can be read from STATE.OUT 601logical :: trdstate=.false. 602! temperature in degrees Kelvin 603real(8) tempk 604 605!--------------------------! 606! mixing variables ! 607!--------------------------! 608! type of mixing to use for the potential 609integer mixtype 610! mixing type description 611character(256) mixdescr 612! adaptive mixing parameters (formerly beta0 and betamax) 613real(8) amixpm(2) 614! subspace dimension for Broyden mixing 615integer mixsdb 616! Broyden mixing parameters alpha and w0 617real(8) broydpm(2) 618 619!----------------------------------------------! 620! charge, moment and current variables ! 621!----------------------------------------------! 622! tolerance for error in total charge 623real(8) epschg 624! total nuclear charge 625real(8) chgzn 626! core charges 627real(8) chgcr(maxspecies) 628! total core charge 629real(8) chgcrtot 630! core leakage charge 631real(8), allocatable :: chgcrlk(:) 632! total valence charge 633real(8) chgval 634! excess charge 635real(8) chgexs 636! total charge 637real(8) chgtot 638! calculated total charge 639real(8) chgcalc 640! interstitial region charge 641real(8) chgir 642! muffin-tin charges 643real(8), allocatable :: chgmt(:) 644! total muffin-tin charge 645real(8) chgmttot 646! effective Wigner radius 647real(8) rwigner 648! total moment 649real(8) momtot(3) 650! total moment magnitude 651real(8) momtotm 652! interstitial region moment 653real(8) momir(3) 654! muffin-tin moments 655real(8), allocatable :: mommt(:,:) 656! total muffin-tin moment 657real(8) mommttot(3) 658! total gauge-invariant current and its magnitude 659real(8) jtot(3),jtotm 660 661!-----------------------------------------! 662! APW and local-orbital variables ! 663!-----------------------------------------! 664! energy step used for numerical calculation of energy derivatives 665real(8) deapwlo 666! maximum allowable APW order 667integer, parameter :: maxapword=4 668! APW order 669integer apword(0:maxlapw,maxspecies) 670! maximum of apword over all angular momenta and species 671integer apwordmax 672! total number of APW coefficients (l, m and order) for each species 673integer lmoapw(maxspecies) 674! polynomial order used for APW radial derivatives 675integer npapw 676! APW initial linearisation energies 677real(8) apwe0(maxapword,0:maxlapw,maxspecies) 678! APW linearisation energies 679real(8), allocatable :: apwe(:,:,:) 680! APW derivative order 681integer apwdm(maxapword,0:maxlapw,maxspecies) 682! apwve is .true. if the linearisation energies are allowed to vary 683logical apwve(maxapword,0:maxlapw,maxspecies) 684! APW radial functions 685real(8), allocatable :: apwfr(:,:,:,:,:) 686! derivate of radial functions at the muffin-tin surface 687real(8), allocatable :: apwdfr(:,:,:) 688! maximum number of local-orbitals 689integer, parameter :: maxlorb=200 690! maximum allowable local-orbital order 691integer, parameter :: maxlorbord=5 692! number of local-orbitals 693integer nlorb(maxspecies) 694! maximum nlorb over all species 695integer nlomax 696! total number of local-orbitals 697integer nlotot 698! local-orbital order 699integer lorbord(maxlorb,maxspecies) 700! maximum lorbord over all species 701integer lorbordmax 702! polynomial order used for local-orbital radial derivatives 703integer nplorb 704! local-orbital angular momentum 705integer lorbl(maxlorb,maxspecies) 706! maximum lorbl over all species 707integer lolmax 708! (lolmax+1)^2 709integer lolmmax 710! local-orbital initial energies 711real(8) lorbe0(maxlorbord,maxlorb,maxspecies) 712! local-orbital energies 713real(8), allocatable :: lorbe(:,:,:) 714! local-orbital derivative order 715integer lorbdm(maxlorbord,maxlorb,maxspecies) 716! lorbve is .true. if the linearisation energies are allowed to vary 717logical lorbve(maxlorbord,maxlorb,maxspecies) 718! local-orbital radial functions 719real(8), allocatable :: lofr(:,:,:,:) 720! band energy search tolerance 721real(8) epsband 722! maximum allowed change in energy during band energy search; enforced only if 723! default energy is less than zero 724real(8) demaxbnd 725! minimum default linearisation energy over all APWs and local-orbitals 726real(8) e0min 727! if autolinengy is .true. then the fixed linearisation energies are set to the 728! Fermi energy minus dlefe 729logical autolinengy 730! difference between linearisation and Fermi energies when autolinengy is .true. 731real(8) dlefe 732! lorbcnd is .true. if conduction state local-orbitals should be added 733logical lorbcnd 734! conduction state local-orbital order 735integer lorbordc 736! excess order of the APW and local-orbital functions 737integer nxoapwlo 738! excess local orbitals 739integer nxlo 740 741!-------------------------------------------! 742! overlap and Hamiltonian variables ! 743!-------------------------------------------! 744! overlap and Hamiltonian matrices sizes at each k-point 745integer, allocatable :: nmat(:,:) 746! maximum nmat over all k-points 747integer nmatmax 748! index to the position of the local-orbitals in the H and O matrices 749integer, allocatable :: idxlo(:,:,:) 750! APW-local-orbital overlap integrals 751real(8), allocatable :: oalo(:,:,:) 752! local-orbital-local-orbital overlap integrals 753real(8), allocatable :: ololo(:,:,:) 754! APW-APW Hamiltonian integrals 755real(8), allocatable :: haa(:,:,:,:,:,:) 756! local-orbital-APW Hamiltonian integrals 757real(8), allocatable :: hloa(:,:,:,:,:) 758! local-orbital-local-orbital Hamiltonian integrals 759real(8), allocatable :: hlolo(:,:,:,:) 760! complex Gaunt coefficient array 761complex(8), allocatable :: gntyry(:,:,:) 762! tefvr is .true. if the first-variational eigenvalue equation is to be solved 763! as a real symmetric problem 764logical tefvr 765! tefvit is .true. if the first-variational eigenvalue equation is to be solved 766! iteratively 767logical tefvit 768! minimum and maximum allowed number of eigenvalue equation iterations 769integer minitefv,maxitefv 770! eigenvalue mixing parameter for iterative solver 771real(8) befvit 772! iterative solver convergence tolerance 773real(8) epsefvit 774 775!--------------------------------------------! 776! eigenvalue and occupancy variables ! 777!--------------------------------------------! 778! number of empty states per atom and spin 779real(8) nempty0 780! number of empty states 781integer nempty 782! number of first-variational states 783integer nstfv 784! number of second-variational states 785integer nstsv 786! smearing type 787integer stype 788! smearing function description 789character(256) sdescr 790! smearing width 791real(8) swidth 792! autoswidth is .true. if the smearing width is to be determined automatically 793logical autoswidth 794! effective mass used in smearing width formula 795real(8) mstar 796! maximum allowed occupancy (1 or 2) 797real(8) occmax 798! convergence tolerance for occupancies 799real(8) epsocc 800! second-variational occupation numbers 801real(8), allocatable :: occsv(:,:) 802! Fermi energy for second-variational states 803real(8) efermi 804! scissor correction applied when computing response functions 805real(8) scissor 806! density of states at the Fermi energy 807real(8) fermidos 808! estimated indirect and direct band gaps 809real(8) bandgap(2) 810! k-points of indirect and direct gaps 811integer ikgap(3) 812! error tolerance for the first-variational eigenvalues 813real(8) evaltol 814! second-variational eigenvalues 815real(8), allocatable :: evalsv(:,:) 816! tevecsv is .true. if second-variational eigenvectors are calculated 817logical tevecsv 818! maximum number of k-point and states indices in user-defined list 819integer, parameter :: maxkst=20 820! number of k-point and states indices in user-defined list 821integer nkstlist 822! user-defined list of k-point and state indices 823integer kstlist(2,maxkst) 824 825!------------------------------! 826! core state variables ! 827!------------------------------! 828! occupancies for core states 829real(8), allocatable :: occcr(:,:) 830! eigenvalues for core states 831real(8), allocatable :: evalcr(:,:) 832! radial wavefunctions for core states 833real(8), allocatable :: rwfcr(:,:,:,:) 834! radial charge density for core states 835real(8), allocatable :: rhocr(:,:,:) 836! spincore is .true. if the core is to be treated as spin-polarised 837logical spincore 838! number of core spin-channels 839integer nspncr 840 841!--------------------------! 842! energy variables ! 843!--------------------------! 844! eigenvalue sum 845real(8) evalsum 846! electron kinetic energy 847real(8) engykn 848! core electron kinetic energy 849real(8) engykncr 850! nuclear-nuclear energy 851real(8) engynn 852! electron-nuclear energy 853real(8) engyen 854! Hartree energy 855real(8) engyhar 856! Coulomb energy (E_nn + E_en + E_H) 857real(8) engycl 858! electronic Coulomb potential energy 859real(8) engyvcl 860! Madelung term 861real(8) engymad 862! exchange-correlation potential energy 863real(8) engyvxc 864! exchange-correlation effective field energy 865real(8) engybxc 866! energy of external global magnetic field 867real(8) engybext 868! exchange energy 869real(8) engyx 870! correlation energy 871real(8) engyc 872! electronic entropy 873real(8) entrpy 874! entropic contribution to free energy 875real(8) engyts 876! total energy 877real(8) engytot 878 879!--------------------------------------------! 880! force, stress and strain variables ! 881!--------------------------------------------! 882! tforce is .true. if force should be calculated 883logical tforce 884! Hellmann-Feynman force on each atom 885real(8), allocatable :: forcehf(:,:) 886! incomplete basis set (IBS) force on each atom 887real(8), allocatable :: forceibs(:,:) 888! total force on each atom 889real(8), allocatable :: forcetot(:,:) 890! previous total force on each atom 891real(8), allocatable :: forcetotp(:,:) 892! maximum force magnitude over all atoms 893real(8) forcemax 894! tfav0 is .true. if the average force should be zero in order to prevent 895! translation of the atomic basis 896logical tfav0 897! average force 898real(8) forceav(3) 899! atomic position optimisation type 900! 0 : no optimisation 901! 1 : unconstrained optimisation 902integer atpopt 903! maximum number of atomic position optimisation steps 904integer maxatpstp 905! default step size parameter for atomic position optimisation 906real(8) tau0atp 907! step size parameters for each atom 908real(8), allocatable :: tauatp(:) 909! number of strain tensors 910integer nstrain 911! current strain tensor 912integer :: istrain=0 913! strain tensors 914real(8) strain(3,3,9) 915! infinitesimal displacement parameter multiplied by the strain tensor for 916! computing the stress tensor 917real(8) deltast 918! symmetry reduced stress tensor components 919real(8) stress(9) 920! previous stress tensor 921real(8) stressp(9) 922! stress tensor component magnitude maximum 923real(8) stressmax 924! lattice vector optimisation type 925! 0 : no optimisation 926! 1 : unconstrained optimisation 927! 2 : iso-volumetric optimisation 928integer latvopt 929! maximum number of lattice vector optimisation steps 930integer maxlatvstp 931! default step size parameter for lattice vector optimisation 932real(8) tau0latv 933! step size for each stress tensor component acting on the lattice vectors 934real(8) taulatv(9) 935 936!-------------------------------! 937! convergence variables ! 938!-------------------------------! 939! maximum number of self-consistent loops 940integer maxscl 941! current self-consistent loop number 942integer iscl 943! Kohn-Sham potential convergence tolerance 944real(8) epspot 945! energy convergence tolerance 946real(8) epsengy 947! force convergence tolerance 948real(8) epsforce 949! stress tensor convergence tolerance 950real(8) epsstress 951 952!----------------------------------------------------------! 953! density of states, optics and response variables ! 954!----------------------------------------------------------! 955! number of energy intervals in the DOS/optics function plot 956integer nwplot 957! fine k-point grid size for integration of functions in the Brillouin zone 958integer ngrkf 959! smoothing level for DOS/optics function plot 960integer nswplot 961! energy interval for DOS/optics function plot 962real(8) wplot(2) 963! maximum angular momentum for the partial DOS plot 964integer lmaxdos 965! dosocc is .true. if the DOS is to be weighted by the occupancy 966logical dosocc 967! dosmsum is .true. if the partial DOS is to be summed over m 968logical dosmsum 969! dosssum is .true. if the partial DOS is to be summed over spin 970logical dosssum 971! number of optical matrix components required 972integer noptcomp 973! required optical matrix components 974integer optcomp(3,27) 975! intraband is .true. if the intraband term is to be added to the optical matrix 976logical intraband 977! lmirep is .true. if the (l,m) band characters should correspond to the 978! irreducible representations of the site symmetries 979logical lmirep 980! spin-quantisation axis in Cartesian coordinates used when plotting the 981! spin-resolved DOS (z-axis by default) 982real(8) sqados(3) 983! q-vector in lattice and Cartesian coordinates for calculating the matrix 984! elements < i,k+q | exp(iq.r) | j,k > 985real(8) vecql(3),vecqc(3) 986! maximum initial-state energy allowed in ELNES transitions 987real(8) emaxelnes 988! structure factor energy window 989real(8) wsfac(2) 990 991!-------------------------------------! 992! 1D/2D/3D plotting variables ! 993!-------------------------------------! 994! number of vertices in 1D plot 995integer nvp1d 996! total number of points in 1D plot 997integer npp1d 998! vertices in lattice coordinates for 1D plot 999real(8), allocatable :: vvlp1d(:,:) 1000! distance to vertices in 1D plot 1001real(8), allocatable :: dvp1d(:) 1002! plot vectors in lattice coordinates for 1D plot 1003real(8), allocatable :: vplp1d(:,:) 1004! distance to points in 1D plot 1005real(8), allocatable :: dpp1d(:) 1006! corner vectors of 2D plot in lattice coordinates 1007real(8) vclp2d(3,0:2) 1008! grid sizes of 2D plot 1009integer np2d(2) 1010! corner vectors of 3D plot in lattice coordinates 1011real(8) vclp3d(3,0:3) 1012! grid sizes of 3D plot 1013integer np3d(3) 1014 1015!----------------------------------------! 1016! OEP and Hartree-Fock variables ! 1017!----------------------------------------! 1018! maximum number of core states over all species 1019integer ncrmax 1020! maximum number of OEP iterations 1021integer maxitoep 1022! OEP step size and fraction of previous potential to be used as starting point 1023real(8) tauoep(2) 1024! magnitude of the OEP residual 1025real(8) resoep 1026! exchange potential and magnetic field 1027real(8), allocatable :: vxmt(:,:),vxir(:) 1028real(8), allocatable :: bxmt(:,:,:),bxir(:,:) 1029! hybrid is .true. if a hybrid functional is to be used 1030logical hybrid,hybrid0 1031! hybrid functional mixing coefficient 1032real(8) hybridc 1033 1034!-------------------------------------------------------------! 1035! response function and perturbation theory variables ! 1036!-------------------------------------------------------------! 1037! |G| cut-off for response functions 1038real(8) gmaxrf 1039! energy cut-off for response functions 1040real(8) emaxrf 1041! number of G-vectors for response functions 1042integer ngrf 1043! number of response function frequencies 1044integer nwrf 1045! complex response function frequencies 1046complex(8), allocatable :: wrf(:) 1047! maximum number of spherical Bessel functions on the coarse radial mesh over 1048! all species 1049integer njcmax 1050 1051!-------------------------------------------------! 1052! Bethe-Salpeter equation (BSE) variables ! 1053!-------------------------------------------------! 1054! number of valence and conduction states for transitions 1055integer nvbse,ncbse 1056! default number of valence and conduction states 1057integer nvbse0,ncbse0 1058! maximum number of extra valence and conduction states 1059integer, parameter :: maxxbse=20 1060! number of extra valence and conduction states 1061integer nvxbse,ncxbse 1062! extra valence and conduction states 1063integer istxbse(maxxbse),jstxbse(maxxbse) 1064! total number of transitions 1065integer nvcbse 1066! size of blocks in BSE Hamiltonian matrix 1067integer nbbse 1068! size of BSE matrix (= 2*nbbse) 1069integer nmbse 1070! index from BSE valence states to second-variational states 1071integer, allocatable :: istbse(:,:) 1072! index from BSE conduction states to second-variational states 1073integer, allocatable :: jstbse(:,:) 1074! index from BSE valence-conduction pair and k-point to location in BSE matrix 1075integer, allocatable :: ijkbse(:,:,:) 1076! BSE Hamiltonian 1077complex(8), allocatable :: hmlbse(:,:) 1078! BSE Hamiltonian eigenvalues 1079real(8), allocatable :: evalbse(:) 1080! if bsefull is .true. then the full BSE Hamiltonian is calculated, otherwise 1081! only the Hermitian block 1082logical bsefull 1083! if hxbse/hdbse is .true. then the exchange/direct term is included in the BSE 1084! Hamiltonian 1085logical hxbse,hdbse 1086 1087!--------------------------! 1088! timing variables ! 1089!--------------------------! 1090! initialisation 1091real(8) timeinit 1092! Hamiltonian and overlap matrix set up 1093real(8) timemat 1094! first-variational calculation 1095real(8) timefv 1096! second-variational calculation 1097real(8) timesv 1098! charge density calculation 1099real(8) timerho 1100! potential calculation 1101real(8) timepot 1102! force calculation 1103real(8) timefor 1104 1105!-----------------------------! 1106! numerical constants ! 1107!-----------------------------! 1108real(8), parameter :: pi=3.1415926535897932385d0 1109real(8), parameter :: twopi=6.2831853071795864769d0 1110real(8), parameter :: fourpi=12.566370614359172954d0 1111! spherical harmonic for l=m=0 1112real(8), parameter :: y00=0.28209479177387814347d0 1113! complex constants 1114complex(8), parameter :: zzero=(0.d0,0.d0) 1115complex(8), parameter :: zone=(1.d0,0.d0) 1116complex(8), parameter :: zi=(0.d0,1.d0) 1117! array of i^l and (-i)^l values 1118integer, private :: l 1119complex(8), parameter :: zil(0:maxlapw)=[((0.d0,1.d0)**l,l=0,maxlapw)] 1120complex(8), parameter :: zilc(0:maxlapw)=[((0.d0,-1.d0)**l,l=0,maxlapw)] 1121! Pauli spin matrices: 1122! sigma_x = ( 0 1 ) sigma_y = ( 0 -i ) sigma_z = ( 1 0 ) 1123! ( 1 0 ) ( i 0 ) ( 0 -1 ) 1124! Planck constant in SI units (exact, CODATA 2018) 1125real(8), parameter :: h_si=6.62607015d-34 1126! reduced Planck constant in SI units 1127real(8), parameter :: hbar_si=h_si/twopi 1128! speed of light in SI units (exact, CODATA 2018) 1129real(8), parameter :: sol_si=299792458d0 1130! speed of light in atomic units (=1/alpha) (CODATA 2018) 1131real(8), parameter :: sol=137.035999084d0 1132! scaled speed of light 1133real(8) solsc 1134! Hartree in SI units (CODATA 2018) 1135real(8), parameter :: ha_si=4.3597447222071d-18 1136! Hartree in eV (CODATA 2018) 1137real(8), parameter :: ha_ev=27.211386245988d0 1138! Hartree in inverse meters 1139real(8), parameter :: ha_im=ha_si/(h_si*sol_si) 1140! Boltzmann constant in SI units (exact, CODATA 2018) 1141real(8), parameter :: kb_si=1.380649d-23 1142! Boltzmann constant in Hartree/kelvin 1143real(8), parameter :: kboltz=kb_si/ha_si 1144! electron charge in SI units (exact, CODATA 2018) 1145real(8), parameter :: e_si=1.602176634d-19 1146! Bohr radius in SI units (CODATA 2018) 1147real(8), parameter :: br_si=0.529177210903d-10 1148! Bohr radius in Angstroms 1149real(8), parameter :: br_ang=br_si*1.d10 1150! atomic unit of magnetic flux density in SI 1151real(8), parameter :: b_si=hbar_si/(e_si*br_si**2) 1152! atomic unit of electric field in SI 1153real(8), parameter :: ef_si=ha_si/(e_si*br_si) 1154! atomic unit of time in SI 1155real(8), parameter :: t_si=hbar_si/ha_si 1156! electron g-factor (CODATA 2018) 1157real(8), parameter :: gfacte=2.00231930436256d0 1158! electron mass in SI (CODATA 2018) 1159real(8), parameter :: em_si=9.1093837015d-31 1160! atomic mass unit in SI (CODATA 2018) 1161real(8), parameter :: amu_si=1.66053906660d-27 1162! atomic mass unit in electron masses 1163real(8), parameter :: amu=amu_si/em_si 1164 1165!---------------------------------! 1166! miscellaneous variables ! 1167!---------------------------------! 1168! code version 1169integer, parameter :: version(3)=[7,2,42] 1170! maximum number of tasks 1171integer, parameter :: maxtasks=40 1172! number of tasks 1173integer ntasks 1174! task array 1175integer tasks(maxtasks) 1176! current task 1177integer task 1178! tlast is .true. if the calculation is on the last self-consistent loop 1179logical tlast 1180! tstop is .true. if the STOP file exists 1181logical tstop 1182! number of self-consistent loops after which STATE.OUT is written 1183integer nwrite 1184! filename extension for files generated by gndstate 1185character(256) :: filext='.OUT' 1186! scratch space path 1187character(256) scrpath 1188! number of note lines 1189integer notelns 1190! notes to include in INFO.OUT 1191character(256), allocatable :: notes(:) 1192 1193end module 1194 1195