~siesta-maint/siesta/trunk : revision 678

5

! or http://www.gnu.org/copyleft/gpl.txt.

6

! See Docs/Contributors.txt for a list of contributors.

7

!

8

subroutine mulliken(iopt,nspin,natoms,nbasistot,maxnh,numh,

8

subroutine mulliken(iopt,natoms,nbasistot,maxnh,numh,

9

. listhptr,listh,s,dm,isa,lasto,iaorb,

10

. iphorb)

11

C ********************************************************************

23

C integer iopt : Work option: 1 = atomic and orbital charges

24

C 2 = 1 + atomic overlap pop.

25

C 3 = 2 + orbital overlap pop.

26

C integer nspin : Number of spin components

27

26

C integer natoms : Number of atoms

28

27

C integer nbasistot : Number of basis orbitals over all nodes

29

28

C integer maxnh : First dimension of d.m. and overlap, and its

32

31

C integer listhptr(nbasis) : Second Control vector of d.m. and overlap

33

32

C integer listh(maxnh) : Third Control vector of d.m. and overlap

34

33

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

35

C real*8 dm(maxnh,nspin) : Density matrix in sparse form

34

C real*8 dm(maxnh,spin%DM) : Density matrix in sparse form

36

35

C integer isa(natoms) : Species index of each atom

37

36

C integer lasto(0:maxa) : Index of last orbital of each atom

38

37

C (lasto(0) = 0)

51

50

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

52

51

use alloc, only : re_alloc, de_alloc

53

52

use siesta_cml

53

use m_spin, only : spin

54

#ifdef MPI

55

use mpi_siesta

56

#endif

57

58

implicit none

59

60

integer

61

. iopt,natoms,nbasistot,maxnh,nspin

60

integer, intent(in) ::

61

. iopt,natoms,nbasistot,maxnh

62

63

integer

64

. numh(*),lasto(0:natoms),listh(maxnh),listhptr(*),

65

. iphorb(*), isa(natoms), iaorb(*)

66

67

real(dp)

68

. dm(maxnh,nspin),s(maxnh)

68

. dm(maxnh,spin%DM),s(maxnh)

69

70

C Internal parameters ..................................................

71

C Number of columns in printout. Must be smaller than 20

84

. config(ncol), itot

85

86

real(dp)

87

. p(ncol),qa,qas(4),qts(4),qtot,qtot_temp,stot,svec(3)

87

. p(ncol),qa,qas(4),qts(4),qtot,qtot_temp,stot,svec(3)

88

real(dp) :: dm_aux(4)

88

89

90

real(dp), pointer :: qo(:)=>null(), qos(:,:)=>null()

90

91

177

178

ii = j - (ib - 1) * ncol

178

179

if (ii .ge. 1 .and. ii .le. imax) then

179

180

p(ii) = 0._dp

180

do is = 1,min(nspin,2)

181

do is = 1,spin%Spinor

181

182

p(ii) = p(ii) + dm(ind,is) * s(ind)

182

183

enddo

183

184

endif

241

242

goto 110

242

243

100 ii = ja - (ib - 1) * ncol

243

244

if (ii .ge. 1 .and. ii .le. imax) then

244

do is = 1,min(nspin,2)

245

do is = 1,spin%Spinor

245

246

p(ii) = p(ii) + dm(ind,is) * s(ind)

246

247

enddo

247

248

endif

271

272

write(6,*)

272

273

write(6,"(a)")'mulliken: Atomic and Orbital Populations:'

273

274

endif

274

if (nspin .le. 2) then

275

if (spin%DM .le. 2) then

275

276

C We only write mulliken to cml file if we are unpolarized or collinear

276

277

C and (currently) don't bother in the non-collinear case.

277

278

if (cml_p) then

279

280

. dictRef='siesta:mulliken',

280

281

. title='Mulliken Population Analysis')

281

282

endif

282

do is = 1,nspin

283

if (nspin .eq. 2) then

283

do is = 1,spin%DM

284

if (spin%DM .eq. 2) then

284

285

if(is .eq. 1) then

285

286

if (ionode) write(6,'(/,a)') 'mulliken: Spin UP '

286

287

if (cml_p) call cmlStartPropertyList(xf=mainXML,

399

400

enddo

400

401

if (cml_p) call cmlEndPropertyList(xf=mainXML) ! mulliken

401

402

elseif (nspin .eq. 4) then

403

call re_alloc( qos, 1, nspin, 1, nbasistot,

403

elseif (spin%DM .ge. 4) then

404

call re_alloc( qos, 1, spin%Grid, 1, nbasistot,

404

405

& 'qos', 'mulliken' )

405

do is = 1,nspin

406

do is = 1,spin%Grid

406

407

qts(is) = 0.0_dp

407

408

do io = 1,nbasistot

408

409

qos(is,io) = 0.0_dp

410

411

enddo

411

412

! For SOC-onsite, dm(:,4) has the wrong sign, but this

412

413

! routine is not called ('moments' is).

413

do is = 1,nspin

414

do io = 1,nbasis

415

call LocalToGlobalOrb(io,Node, Nodes, itot)

416

do in = 1,numh(io)

414

do io = 1,nbasis

415

call LocalToGlobalOrb(io,Node, Nodes, itot)

416

do in = 1,numh(io)

417

ind = listhptr(io)+in

418

jo = listh(ind)

419

qos(is,itot) = qos(is,itot) + dm(ind,is)*s(ind)/2

420

if ( jo <= nbasistot )

421

. qos(is,jo) = qos(is,jo) + dm(ind,is)*s(ind)/2

422

enddo

423

enddo

419

dm_aux(1:4) = dm(ind,1:4)

420

if (spin%SO) then

421

dm_aux(3) = 0.5_dp*(dm_aux(3)+dm(ind,7))

422

dm_aux(4) = 0.5_dp*(dm_aux(4)+dm(ind,8))

423

endif

424

do is = 1,spin%Grid

425

qos(is,itot) = qos(is,itot) + dm_aux(is)*s(ind)/2

426

if ( jo <= nbasistot )

427

. qos(is,jo) = qos(is,jo) + dm_aux(is)*s(ind)/2

428

enddo

429

enddo

424

430

enddo

425

431

#ifdef MPI

426

call re_alloc( qb, 1, nspin*nbasistot,

432

call re_alloc( qb, 1, spin%Grid*nbasistot,

427

433

& 'qb', 'mulliken' )

428

call MPI_Reduce(qos(1,1), qb(1), nspin*nbasistot,

434

call MPI_Reduce(qos(1,1), qb(1), spin%Grid*nbasistot,

429

435

& Mpi_Double_Precision, MPI_Sum,

430

436

& 0, MPI_Comm_world, MPIerror)

431

437

itot = 0

432

438

do io = 1 , nbasistot

433

do is = 1 , nspin

439

do is = 1 , spin%Grid

434

440

itot = itot + 1

435

441

qos(is, io) = qb(itot)

436

442

end do

462

468

if ( isa(ia) /= ispec ) cycle

463

469

464

470

! Initialize atomic charges (per spin)

465

qas(1:nspin) = 0._dp

471

qas(1:spin%Grid) = 0._dp

466

472

do io = lasto(ia-1)+1, lasto(ia)

467

473

468

474

C DSP, Writing symmetries for each orbital.

471

477

config(1) = cnfigfio(ispec,iphorb(io))

472

478

473

479

! Reduce total charge per-atom

474

do is = 1 , nspin

480

do is = 1 , spin%Grid

475

481

qas(is) = qas(is) + qos(is,io)

476

482

enddo

477

483

478

call spnvec( nspin, qos(1:nspin,io),

484

call spnvec( spin%Grid, qos(1:spin%Grid,io),

479

485

$ qtot, stot, svec)

480

486

if ( IONode ) then

481

487

write(6,'(i5,i3,tr1,i1,a7,2f10.5,3x,3f8.3)')

485

491

486

492

end do ! atomic orbitals

487

493

488

call spnvec( nspin, qas, qtot, stot, svec )

494

call spnvec( spin%Grid, qas, qtot, stot, svec )

489

495

if ( IONode ) then

490

496

write(6,'(i5,5x,a5,2x,2f10.5,3x,3f8.3,/)')

491

497

. ia, 'Total', qtot, stot, svec

492

498

end if

493

499

494

500

! Reduce total charge in system

495

do is = 1 , nspin

501

do is = 1 , spin%Grid

496

502

qts(is) = qts(is) + qas(is)

497

503

end do

498

504

499

505

end do ! loop atoms

500

506

501

call spnvec( nspin, qts, qtot, stot, svec )

507

call spnvec( spin%Grid, qts, qtot, stot, svec )

502

508

if ( IONode ) then

503

509

write(6,'(a,/,a5,12x,2f10.5,3x,3f8.3,/)')

504

510

. repeat('-',64), 'Total', qtot, stot, svec

508

514

509

515

call de_alloc( qos, 'qos', 'mulliken' )

510

516

511

end if ! nspin .eq. 4

517

end if ! spin%DM .ge. 4

512

518

end if

513

519

514

520

C ...................

521

527

11 format(i12,20(1x,f7.3))

522

528

15 format(i4,1x,i1,a7,20(1x,f8.3))

523

529

return

524

end

530

end subroutine mulliken

525

531

526

532

527

533

582

588

write(6,*) 'spnvec: ERROR: invalid argument ns =', ns

583

589

return

584

590

endif

585

end

591

end subroutine spnvec