~siesta-pseudos-bases/siesta/trunk-psml

« back to all changes in this revision

Viewing changes to Src/bands.F

• browse files at revision 605
• compare with another revision
• download diff
• download tarball
• view history from revision 605

Show diffs side-by-side

added

removed

Src/bands.F

296

296

      end subroutine initbands

297

297

 

298

298

 

299

 

      subroutine bands( no_s, nspin, maxspn, no_u, no_l, maxnh, maxk,

 

299

      subroutine bands( no_s, spin,

 

300

     .                  no_u, no_l,maxnh, maxk,

300

301

     .                  numh, listhptr, listh, H, S, ef, xij, indxuo,

301

302

     .                  writeb, nk, kpoint, ek, occtol, getPSI )

302

303

 

305

306

C Written by J.Soler, August 1997 and August 1998.

306

307

C Initialisation moved into a separate routine, JDG Jan 2000.

307

308

C WFS options by A. Garcia, April 2012

 

309

C Cleaned for tSpin by N. Papior, August 2019

308

310

 

309

311

C **************************** INPUT **********************************

310

312

C integer no_s                : Number of basis orbitals in supercell

311

 

C integer nspin               : Number of spin components

312

 

C integer maxspn              : Maximum number of spin components

313

 

C                               nspin /= maxspn for spin-orbit

 

313

C type(tSpin) spin            : Containing all spin-information

314

314

C integer maxnh               : Maximum number of orbitals interacting

315

315

C                               with any orbital

316

316

C integer maxk                : Last dimension of kpoint and ek

317

 

C integer numh(no_l)           : Number of nonzero elements of each row

 

317

C integer numh(no_l)          : Number of nonzero elements of each row

318

318

C                               of hamiltonian matrix

319

 

C integer listhptr(no_l)       : Pointer to start of each row of the

 

319

C integer listhptr(no_l)      : Pointer to start of each row of the

320

320

C                               hamiltonian matrix

321

321

C integer listh(maxlh)        : Nonzero hamiltonian-matrix element

322

322

C                               column indexes for each matrix row

323

 

C real*8  H(maxnh,nspin)      : Hamiltonian in sparse form

 

323

C real*8  H(maxnh,spin%H)     : Hamiltonian in sparse form

324

324

C real*8  S(maxnh)            : Overlap in sparse form

325

325

C real*8  ef                  : Fermi energy

326

326

C real*8  xij(3,maxnh)        : Vectors between orbital centers (sparse)

327

327

C                               (not used if only gamma point)

328

328

C integer no_u                : First dimension of ek

329

 

C integer no_l               : Second dimension of H and S

330

 

C integer indxuo(no_s)          : Index of equivalent orbital in unit cell

 

329

C integer no_l                : Second dimension of H and S

 

330

C integer indxuo(no_s)        : Index of equivalent orbital in unit cell

331

331

C                               Unit cell orbitals must be the first in

332

332

C                               orbital lists, i.e. indxuo.le.no_l, with

333

333

C                               no_l the number of orbitals in unit cell

334

334

C real*8  ef                  : Fermi energy

335

335

C logical writeb              : This routine must write bands?

336

 

C integer no_u              : Total number of orbitals in unit cell

 

336

C integer no_u                : Total number of orbitals in unit cell

337

337

C integer nk                  : Number of band k points

338

338

C real*8  kpoint(3,maxk)      : k point vectors

339

339

C real*8  occtol              : Occupancy threshold for DM build

340

340

C *************************** OUTPUT **********************************

341

 

C real*8  ek(no_u,maxspn,maxk) : Eigenvalues

 

341

C real*8  ek(no_u,spin%spinor,maxk) : Eigenvalues

342

342

C *************************** UNITS ***********************************

343

343

C Lengths in atomic units (Bohr).

344

344

C k vectors in reciprocal atomic units.

366

366

      use atmfuncs,     only : symfio, cnfigfio, labelfis, nofis

367

367

      use writewave,    only : wfs_filename

368

368

 

369

 

      use m_spin,       only: spin 

 

369

      use t_spin, only: tSpin

370

370

 

371

371

      use m_diag_option, only: ParallelOverK, Serial

372

372

      use m_diag, only: diag_init

374

374

 

375

375

      implicit          none

376

376

 

377

 

      integer  ::   nspin

378

 

      integer  ::   maxk, maxnh, maxspn, no_u, no_l, nk, no_s,

 

377

      type(tSpin), intent(in) :: spin

 

378

      integer  ::   maxk, maxnh, no_u, no_l, nk, no_s,

379

379

     .              indxuo(no_s), listh(maxnh), numh(no_l),

380

 

     .              listhptr(*), nuo

 

380

     .              listhptr(*)

381

381

      logical  ::   writeb

382

 

      real(dp) ::   ef, ek(no_u,maxspn,maxk),

383

 

     .              H(maxnh,nspin), kpoint(3,maxk),

 

382

      real(dp) ::   ef, ek(no_u,spin%spinor,maxk),

 

383

     .              H(maxnh,spin%H), kpoint(3,maxk),

384

384

     .              S(maxnh), xij(3,maxnh), occtol

385

385

      logical, intent(in) :: getPSI

386

386

 

394

394

      logical, parameter :: fixspin = .false.

395

395

 

396

396

      integer :: ik, il, io, ispin, iu, iu_wfs, iuo, naux, nhs, j

397

 

 

398

397

      logical :: SaveParallelOverK

399

398

      

400

399

      real(dp)

405

404

      real(dp), dimension(:), pointer :: aux => null()

406

405

 

407

406

      parameter ( eV = 1.d0 / 13.60580d0 )

 

407

 

408

408

      save Dnew, Enew, e1, e2, qk, qtot, temp, wk

409

 

 

410

409

      data Dnew, Enew, e1, e2, qk, qtot, temp, wk /8*0.d0/

411

410

 

412

 

C Get number of local orbitals

413

 

#ifdef MPI

414

 

      call GetNodeOrbs(no_u,Node,Nodes,nuo)

415

 

#else

416

 

      nuo = no_u

417

 

#endif

418

 

 

419

411

C Start time counter

420

412

      call timer( 'bands', 1 )

421

413

 

430

422

      endif

431

423

 

432

424

C     Allocate local arrays - only aux is relevant here

433

 

      if ( nspin >=4 ) then

434

 

         nhs  = 2 * (2*no_u) * (2*no_l)

 

425

 

 

426

      if ( spin%Grid == 4 ) then

 

427

         nhs = 2 * (2*no_u) * (2*no_l)

435

428

      else

436

429

         nhs = 2 * no_u*no_l

437

430

      endif

438

 

      naux  = 2*no_u*5

 

431

      naux = 2*no_u*5

439

432

      call allocDenseMatrix(nhs, nhs, nhs)

440

433

      call re_alloc( aux, 1, naux, 'aux', 'bands' )

441

434

 

448

441

        rewind (iu_wfs)

449

442

 

450

443

        write(iu_wfs) nk, .false. ! nk, Gamma, same file-format in WFS as for Gamma-point

451

 

        write(iu_wfs) nspin

 

444

        write(iu_wfs) spin%Grid

452

445

        write(iu_wfs) no_u

453

446

        write(iu_wfs) (iaorb(j),labelfis(isa(iaorb(j))),

454

447

     .            iphorb(j), cnfigfio(isa(iaorb(j)),iphorb(j)),

458

451

      endif

459

452

 

460

453

C Find the band energies

461

 

      if ( spin%none .or. spin%Col ) then

 

454

      if ( spin%Grid < 4 ) then

462

455

C fixspin and qs are not used in diagk, since getD=.false. ...

463

456

        qs(1) = 0.0_dp

464

457

        qs(2) = 0.0_dp

475

468

           call diag_init()

476

469

        end if

477

470

 

478

 

        call diagk( nspin, no_l, no_s, maxspn, 

 

471

        call diagk( spin%H, no_l, no_s, spin%spinor, 

479

472

     .              maxnh, maxnh, 

480

473

     .              no_u, numh, listhptr, listh, numh, listhptr,

481

474

     .              listh, H, S, getD, getPSI, fixspin, qtot, qs, temp,

488

481

        if ( ParallelOverK ) Serial = .true.

489

482

 

490

483

      elseif ( spin%NCol ) then

491

 

         if (getPSI) then

492

 

            if (node==0) then

493

 

               write(6,*) "No WFS for non-colinear spin, yet..."

494

 

            endif

495

 

            RETURN

496

 

         endif

497

484

        call diag2k(no_l, no_s, maxnh, maxnh, no_u,

498

485

     .              numh, listhptr, listh, numh, listhptr,

499

 

     .              listh, H, S, getD, qtot, temp, e1, e2,

 

486

     .              listh, H, S, getD, getPSI, qtot, temp, e1, e2,

500

487

     .              xij, indxuo, nk, kpoint, wk,

501

488

     .              ek, qk, Dnew, Enew, ef, Entropy,

502

489

     .              Haux, Saux, psi, Haux, Saux, aux,

503

490

     .              no_u, occtol, 1, no_u )

504

491

 

505

492

      elseif ( spin%SO ) then

506

 

        if (getPSI) then

507

 

           if (node==0) then

508

 

              write(6,*) "No WFS for spin-orbit, yet..."

509

 

           endif

510

 

           return

511

 

        end if

512

 

        call diag3k(nuo, no_s, maxnh, maxnh, no_u, numh,

 

493

        call diag3k(no_l, no_s, maxnh, maxnh, no_u, numh,

513

494

     .              listhptr, listh, numh, listhptr, listh,

514

495

     .              H, S, getD, getPSI, qtot, temp, e1, e2, xij,

515

496

     .              indxuo, nk, kpoint, wk, ek, qk, Dnew,

517

498

     .              Haux, Saux, psi, Haux, Saux, aux,

518

499

     .              no_u, occtol, 1, no_u )

519

500

      else

520

 

        call die( 'bands: ERROR: incorrect value of nspin')

 

501

        call die( 'bands: ERROR: incorrect value of spin%H')

521

502

      endif

522

503

 

523

504

C Write bands

540

521

     .      path = path + sqrt( (kpoint(1,ik)-kpoint(1,ik-1))**2 +

541

522

     .                          (kpoint(2,ik)-kpoint(2,ik-1))**2 +

542

523

     .                          (kpoint(3,ik)-kpoint(3,ik-1))**2 )

543

 

          do ispin = 1,maxspn

 

524

          do ispin = 1, spin%spinor

544

525

            do io = 1, no_u

545

526

              emax = max( emax, ek(io,ispin,ik) )

546

527

              emin = min( emin, ek(io,ispin,ik) )

552

533

        write(iu,*) emin/eV, emax/eV

553

534

 

554

535

C Write eigenvalues

555

 

        if ( nspin <= 2 ) then

556

 

          write(iu,*) no_u, maxspn, nk

557

 

        else if ( nspin <= 8 ) then

558

 

          write(iu,*) 2*no_u, 1, nk

559

 

        else

560

 

          call die( 'bands: ERROR: incorrect value of nspin' )

561

 

        end if

 

536

        if ( spin%Grid == 4 ) then

 

537

          write(iu,*) 2*no_u, 1, nk     ! A single spin channel, double number of bands

 

538

        else 

 

539

          write(iu,*) no_u, spin%spinor, nk

 

540

        endif

 

541

 

 

542

        ! For non-collinear calculations, the ek array will contain the first no_u

 

543

        ! bands in ek(:,1,ik), and the next no_u bands in ek(:,2,ik)   

562

544

        path = 0.d0

563

545

        do ik = 1,nk

564

546

          if (nlines .ge. 0) then

568

550

     .                            (kpoint(3,ik)-kpoint(3,ik-1))**2 )

569

551

            write(iu,'(f10.6,10f12.4,/,(10x,10f12.4))')

570

552

     .         path, ((ek(io,ispin,ik)/eV,io=1,no_u),

571

 

     .         ispin=1,maxspn)

 

553

     .         ispin=1,spin%spinor)

572

554

          else

573

555

            write(iu,'(3f9.5,8f12.4,/,(27x,8f12.4))')

574

556

     .        kpoint(1,ik), kpoint(2,ik), kpoint(3,ik),

575

 

     .        ((ek(io,ispin,ik)/eV,io=1,no_u),ispin=1,maxspn)

 

557

     .        ((ek(io,ispin,ik)/eV,io=1,no_u),ispin=1,spin%spinor)

576

558

          endif

577

559

        enddo

578

560

 

• Older »

Loggerhead is a web-based interface for Breezy
Version: 2.0.1