~nickpapior/siesta/trunk-buds-format0.92 : revision 432.2.109

1

<?xml version="1.0" encoding="utf-8"?>

2

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"

3

"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">

4

5

<head>

6

<title>Literate version of PEXSI LDOS calculator</title>

7

8

9

10

11

12

<![CDATA[/*><!--*/

13

.title { text-align: center; }

14

.todo { font-family: monospace; color: red; }

15

.done { color: green; }

16

.tag { background-color: #eee; font-family: monospace;

17

padding: 2px; font-size: 80%; font-weight: normal; }

18

.timestamp { color: #bebebe; }

19

.timestamp-kwd { color: #5f9ea0; }

20

.right { margin-left: auto; margin-right: 0px; text-align: right; }

21

.left { margin-left: 0px; margin-right: auto; text-align: left; }

22

.center { margin-left: auto; margin-right: auto; text-align: center; }

23

.underline { text-decoration: underline; }

24

#postamble p, #preamble p { font-size: 90%; margin: .2em; }

25

p.verse { margin-left: 3%; }

26

pre {

27

border: 1px solid #ccc;

28

box-shadow: 3px 3px 3px #eee;

29

padding: 8pt;

30

font-family: monospace;

31

overflow: auto;

32

margin: 1.2em;

33

}

34

pre.src {

35

position: relative;

36

overflow: visible;

37

padding-top: 1.2em;

38

}

39

pre.src:before {

40

display: none;

41

position: absolute;

42

background-color: white;

43

top: -10px;

44

right: 10px;

45

padding: 3px;

46

border: 1px solid black;

47

}

48

pre.src:hover:before { display: inline;}

49

pre.src-sh:before { content: 'sh'; }

50

pre.src-bash:before { content: 'sh'; }

51

pre.src-emacs-lisp:before { content: 'Emacs Lisp'; }

52

pre.src-R:before { content: 'R'; }

53

pre.src-perl:before { content: 'Perl'; }

54

pre.src-java:before { content: 'Java'; }

55

pre.src-sql:before { content: 'SQL'; }

56

57

table { border-collapse:collapse; }

58

caption.t-above { caption-side: top; }

59

caption.t-bottom { caption-side: bottom; }

60

td, th { vertical-align:top; }

61

th.right { text-align: center; }

62

th.left { text-align: center; }

63

th.center { text-align: center; }

64

td.right { text-align: right; }

65

td.left { text-align: left; }

66

td.center { text-align: center; }

67

dt { font-weight: bold; }

68

.footpara:nth-child(2) { display: inline; }

69

.footpara { display: block; }

70

.footdef { margin-bottom: 1em; }

71

.figure { padding: 1em; }

72

.figure p { text-align: center; }

73

.inlinetask {

74

padding: 10px;

75

border: 2px solid gray;

76

margin: 10px;

77

background: #ffffcc;

78

}

79

#org-div-home-and-up

80

{ text-align: right; font-size: 70%; white-space: nowrap; }

81

textarea { overflow-x: auto; }

82

.linenr { font-size: smaller }

83

.code-highlighted { background-color: #ffff00; }

84

.org-info-js_info-navigation { border-style: none; }

85

#org-info-js_console-label

86

{ font-size: 10px; font-weight: bold; white-space: nowrap; }

87

.org-info-js_search-highlight

88

{ background-color: #ffff00; color: #000000; font-weight: bold; }

89

/*]]>*/-->

90

</style>

91

92

/*

93

@licstart The following is the entire license notice for the

94

JavaScript code in this tag.

95

96

97

98

The JavaScript code in this tag is free software: you can

99

redistribute it and/or modify it under the terms of the GNU

100

General Public License (GNU GPL) as published by the Free Software

101

Foundation, either version 3 of the License, or (at your option)

102

any later version. The code is distributed WITHOUT ANY WARRANTY;

103

without even the implied warranty of MERCHANTABILITY or FITNESS

104

FOR A PARTICULAR PURPOSE. See the GNU GPL for more details.

105

106

As additional permission under GNU GPL version 3 section 7, you

107

may distribute non-source (e.g., minimized or compacted) forms of

108

that code without the copy of the GNU GPL normally required by

109

section 4, provided you include this license notice and a URL

110

through which recipients can access the Corresponding Source.

111

112

113

@licend The above is the entire license notice

114

for the JavaScript code in this tag.

115

*/

116

<![CDATA[/*><!--*/

117

function CodeHighlightOn(elem, id)

118

{

119

var target = document.getElementById(id);

120

if(null != target) {

121

elem.cacheClassElem = elem.className;

122

elem.cacheClassTarget = target.className;

123

target.className = "code-highlighted";

124

elem.className = "code-highlighted";

125

}

126

}

127

function CodeHighlightOff(elem, id)

128

{

129

var target = document.getElementById(id);

130

if(elem.cacheClassElem)

131

elem.className = elem.cacheClassElem;

132

if(elem.cacheClassTarget)

133

target.className = elem.cacheClassTarget;

134

}

135

/*]]>*///-->

136

</script>

137

</head>

138

<body>

139

140

<h1 class="title">Literate version of PEXSI LDOS calculator</h1>

141

142

<h2>Table of Contents</h2>

143

144

<ul>

145

<li><a href="#sec-1">1. Introduction</a></li>

146

<li><a href="#sec-2">2. Parent routine with dispatch to Siesta post-processing</a></li>

147

<li><a href="#sec-3">3. PEXSI LDOS routine</a>

148

<ul>

149

<li><a href="#sec-3-1">3.1. Main structure</a></li>

150

<li><a href="#sec-3-2">3.2. Routine header</a>

151

<ul>

152

<li><a href="#sec-3-2-1">3.2.1. Used modules</a></li>

153

</ul>

154

</li>

155

<li><a href="#sec-3-3">3.3. Routine variables</a></li>

156

<li><a href="#sec-3-4">3.4. Define communicators</a></li>

157

<li><a href="#sec-3-5">3.5. Re-distribute matrices</a></li>

158

<li><a href="#sec-3-6">3.6. Set options</a></li>

159

<li><a href="#sec-3-7">3.7. Prepare plan</a></li>

160

<li><a href="#sec-3-8">3.8. Load H and S matrices</a></li>

161

<li><a href="#sec-3-9">3.9. Factorization</a></li>

162

<li><a href="#sec-3-10">3.10. Compute the LDOS</a></li>

163

<li><a href="#sec-3-11">3.11. Copy information to Siesta side</a></li>

164

<li><a href="#sec-3-12">3.12. Clean up</a></li>

165

<li><a href="#sec-3-13">3.13. Support routines</a>

166

<ul>

167

<li><a href="#sec-3-13-1">3.13.1. Row and column partition of npPerPole</a></li>

168

<li><a href="#sec-3-13-2">3.13.2. Error dispatcher</a></li>

169

</ul>

170

</li>

171

</ul>

172

</li>

173

</ul>

174

</div>

175

</div>

176

177

-*- mode: Org -*-

178

179

180

181

<h2 id="sec-1">1 Introduction</h2>

182

183

184

This is a spin-polarized version of the SIESTA-PEXSI interface for the calculation of

185

the LDOS (which is a "partial" DM computed in a restricted energy window).

186

Module structure:

187

188

189

190

191

<pre class="src src-f90">MODULE m_pexsi_local_dos

192

private

193

public :: pexsi_local_dos

194

195

CONTAINS

196

<<siesta-side-parent-routine>>

197

<<pexsi-ldos-routine>>

198

end module m_pexsi_local_dos

199

</pre>

200

</div>

201

</div>

202

</div>

203

204

205

<h2 id="sec-2">2 Parent routine with dispatch to Siesta post-processing</h2>

206

207

208

This routine defines the energy window (<code>energy</code> and <code>broadening</code>),

209

and calls the PEXSI routine to compute the partial DM, which is then

210

passed to the <code>dhscf</code> machinery for the calculation of the

211

corresponding charge, which is just the LDOS needed. The output file

212

will have extension <code>.LDSI</code>.

213

214

215

216

217

<pre class="src src-f90">subroutine pexsi_local_dos( )

218

use m_energies

219

220

use sparse_matrices

221

USE siesta_options

222

use siesta_geom

223

use atomlist, only: indxuo, indxua

224

use atomlist, only: qtot, no_u, no_l

225

use atomlist, only: iphorb

226

use atomlist, only: datm, no_s, iaorb

227

use fdf

228

use files, only : slabel

229

use files, only : filesOut_t ! type for output file names

230

use parallel, only: SIESTA_worker

231

use m_ntm

232

use m_forces, only: fa

233

use m_spin, only: nspin, qs, efs

234

use m_dhscf, only: dhscf

235

236

implicit none

237

238

integer :: dummy_iscf = 1

239

real(dp) :: dummy_str(3,3), dummy_strl(3,3) ! for dhscf call

240

real(dp) :: dummy_dipol(3)

241

242

real(dp) :: factor, g2max

243

real(dp) :: energy, broadening

244

245

type(filesOut_t) :: filesOut ! blank output file names

246

247

energy = fdf_get('PEXSI.LocalDOS.Energy',0.0_dp,"Ry")

248

broadening = fdf_get('PEXSI.LocalDOS.Broadening',0.01_dp,"Ry")

249

250

! Note that we re-use Dscf, so it will be obliterated

251

call get_LDOS_SI( no_u, no_l, nspin, &

252

maxnh, numh, listh, H, S, &

253

Dscf, energy, broadening)

254

255

if (SIESTA_worker) then

256

!Find the LDOS in the real space mesh

257

filesOut%rho = trim(slabel) // '.LDSI'

258

g2max = g2cut

259

260

! There is too much clutter here, because dhscf is

261

! a "kitchen-sink" routine that does too many things

262

263

call dhscf( nspin, no_s, iaorb, iphorb, no_l, &

264

no_u, na_u, na_s, isa, xa_last, indxua, &

265

ntm, 0, 0, 0, filesOut, &

266

maxnh, numh, listhptr, listh, Dscf, Datm, maxnh, H, &

267

Enaatm, Enascf, Uatm, Uscf, DUscf, DUext, Exc, Dxc, &

268

dummy_dipol, dummy_str, fa, dummy_strl )

269

endif

270

271

END subroutine pexsi_local_dos

272

</pre>

273

</div>

274

</div>

275

</div>

276

277

278

279

<h2 id="sec-3">3 PEXSI LDOS routine</h2>

280

281

</div><div id="outline-container-sec-3-1" class="outline-3">

282

<h3 id="sec-3-1">3.1 Main structure</h3>

283

284

285

This is the main structure of the LDOS routine.

286

287

288

289

290

291

<<routine-variables>>

292

! -------- for serial compilation

293

#ifndef MPI

294

call die("PEXSI LDOS needs MPI")

295

#else

296

<<define-communicators>>

297

<<re-distribute-matrices>>

298

<<set-options>>

299

<<prepare-plan>>

300

<<load-hs-matrices>>

301

<<factorization>>

302

<<compute-ldos>>

303

<<copy-to-siesta-side>>

304

<<clean-up>>

305

#endif

306

307

CONTAINS

308

309

<<support-routines>>

310

311

end subroutine get_LDOS_SI

312

</pre>

313

</div>

314

</div>

315

</div>

316

317

318

319

<h3 id="sec-3-2">3.2 Routine header</h3>

320

321

322

323

<pre class="src src-f90">subroutine get_LDOS_SI( no_u, no_l, nspin_in, &

324

maxnh, numh, listh, H, S, &

325

LocalDOSDM, energy, broadening)

326

327

<<used-modules>>

328

329

implicit none

330

331

integer, intent(in) :: maxnh, no_u, no_l, nspin_in

332

integer, intent(in), target :: listh(maxnh), numh(no_l)

333

real(dp), intent(in), target :: H(maxnh,nspin_in), S(maxnh)

334

real(dp), intent(in) :: energy, broadening

335

real(dp), intent(out) :: LocalDOSDM(maxnh,nspin_in)

336

</pre>

337

</div>

338

</div>

339

340

341

<h4 id="sec-3-2-1">3.2.1 Used modules</h4>

342

343

344

345

<pre class="src src-f90"> use precision, only : dp

346

use fdf

347

use units, only: eV, pi

348

use sys, only: die

349

use m_mpi_utils, only: broadcast

350

use parallel, only : SIESTA_worker, BlockSize

351

use parallel, only : SIESTA_Group, SIESTA_Comm

352

use m_redist_spmatrix, only: aux_matrix, redistribute_spmatrix

353

use class_Distribution

354

use alloc, only: re_alloc, de_alloc

355

#ifdef MPI

356

use mpi_siesta

357

#endif

358

use iso_c_binding

359

use f_ppexsi_interface, only: f_ppexsi_options

360

use f_ppexsi_interface, only: f_ppexsi_plan_finalize

361

use f_ppexsi_interface, only: f_ppexsi_plan_initialize

362

use f_ppexsi_interface, only: f_ppexsi_selinv_complex_symmetric_matrix

363

use f_ppexsi_interface, only: f_ppexsi_load_real_symmetric_hs_matrix

364

use f_ppexsi_interface, only: f_ppexsi_set_default_options

365

use f_ppexsi_interface, &

366

only: f_ppexsi_symbolic_factorize_complex_symmetric_matrix

367

</pre>

368

</div>

369

</div>

370

</div>

371

</div>

372

373

374

<h3 id="sec-3-3">3.3 Routine variables</h3>

375

376

377

The local variables for the routine must be declared in a certain

378

place for the compiler, but it is more clear to introduce them as they

379

are needed. The <code>routine-variables</code> noweb-ref will be used for this

380

throughout this document.

381

382

383

384

385

<pre class="src src-f90">integer :: ih, i

386

integer :: info

387

logical :: write_ok

388

!------------

389

external :: timer

390

</pre>

391

</div>

392

</div>

393

</div>

394

395

396

<h3 id="sec-3-4">3.4 Define communicators</h3>

397

398

399

<code>World_Comm</code>, which is in principle set to Siesta's copy of

400

<code>MPI_Comm_World</code>, is the global communicator.

401

402

403

404

405

<pre class="src src-f90">integer :: World_Comm, mpirank, ierr

406

!

407

integer :: norbs

408

integer :: nspin

409

</pre>

410

</div>

411

412

413

414

<pre class="src src-f90">!

415

! Our global communicator is a duplicate of the passed communicator

416

!

417

call MPI_Comm_Dup(true_MPI_Comm_World, World_Comm, ierr)

418

call mpi_comm_rank( World_Comm, mpirank, ierr )

419

420

call timer("pexsi-ldos", 1)

421

422

if (SIESTA_worker) then

423

! rename some intent(in) variables, which are only

424

! defined for the Siesta subset

425

norbs = no_u

426

nspin = nspin_in

427

endif

428

!

429

call broadcast(norbs,comm=World_Comm)

430

call broadcast(nspin,World_Comm)

431

</pre>

432

</div>

433

434

435

Now we need to define the Siesta distribution object and the

436

communicator and distribution object for the first team of PEXSI

437

workers, for the purposes of re-distribution of the relevant

438

matrices.

439

440

441

442

For spin, things are a bit more complicated. We need to make sure that

443

the distributions are defined and known to all processors (via actual

444

ranks) with respect to the same reference bridge communicator. For

445

now, this is WorldComm.

446

447

448

449

450

<pre class="src src-f90">integer :: PEXSI_Pole_Group, PEXSI_Spatial_Group, World_Group

451

integer, allocatable :: pexsi_pole_ranks_in_world(:)

452

integer :: nworkers_SIESTA

453

integer, allocatable :: siesta_ranks_in_world(:)

454

integer :: PEXSI_Pole_Group_in_World

455

integer, allocatable :: PEXSI_Pole_ranks_in_World_Spin(:,:)

456

integer :: PEXSI_Pole_Comm, PEXSI_Spatial_Comm, PEXSI_Spin_Comm

457

integer :: numNodesTotal

458

integer :: npPerPole

459

logical :: PEXSI_worker

460

!

461

type(Distribution) :: dist1

462

type(Distribution), allocatable, target :: dist2_spin(:)

463

type(Distribution), pointer :: dist2

464

integer :: pbs, color, spatial_rank, spin_rank

465

</pre>

466

</div>

467

468

469

Define the Siesta distribution. Note that this is known to all nodes.

470

471

472

473

474

<pre class="src src-f90">call MPI_Comm_Group(World_Comm,World_Group, ierr)

475

call MPI_Group_Size(SIESTA_Group, nworkers_SIESTA, ierr)

476

allocate(siesta_ranks_in_world(nworkers_SIESTA))

477

call MPI_Group_translate_ranks( SIESTA_Group, nworkers_SIESTA, &

478

(/ (i,i=0,nworkers_SIESTA-1) /), &

479

World_Group, siesta_ranks_in_world, ierr )

480

call newDistribution(dist1,World_Comm,siesta_ranks_in_world, &

481

TYPE_BLOCK_CYCLIC,BlockSize,"bc dist")

482

deallocate(siesta_ranks_in_world)

483

call MPI_Barrier(World_Comm,ierr)

484

</pre>

485

</div>

486

487

488

For possibly spin-polarized calculations, we split the communicator.

489

Note that we give the user the option of requesting more processors

490

per pole for the LDOS calculation.

491

492

493

494

495

<pre class="src src-f90">call mpi_comm_size( World_Comm, numNodesTotal, ierr )

496

497

npPerPole = fdf_get("PEXSI.np-per-pole",4)

498

npPerPole = fdf_get("PEXSI.LocalDOS.np-per-pole",npPerPole)

499

if (nspin*npPerPole > numNodesTotal) &

500

call die("PEXSI.np-per-pole is too big for MPI size")

501

502

! "Row" communicator for independent PEXSI operations on each spin

503

! The name refers to "spatial" degrees of freedom.

504

color = mod(mpirank,nspin) ! {0,1} for nspin = 2, or {0} for nspin = 1

505

call MPI_Comm_Split(World_Comm, color, mpirank, PEXSI_Spatial_Comm, ierr)

506

507

! "Column" communicator for spin reductions

508

color = mpirank/nspin

509

call MPI_Comm_Split(World_Comm, color, mpirank, PEXSI_Spin_Comm, ierr)

510

511

! Group and Communicator for first-pole team of PEXSI workers

512

!

513

call MPI_Comm_Group(PEXSI_Spatial_Comm, PEXSI_Spatial_Group, Ierr)

514

call MPI_Group_incl(PEXSI_Spatial_Group, npPerPole, &

515

(/ (i,i=0,npPerPole-1) /),&

516

PEXSI_Pole_Group, Ierr)

517

call MPI_Comm_create(PEXSI_Spatial_Comm, PEXSI_Pole_Group,&

518

PEXSI_Pole_Comm, Ierr)

519

520

521

call mpi_comm_rank( PEXSI_Spatial_Comm, spatial_rank, ierr )

522

call mpi_comm_rank( PEXSI_Spin_Comm, spin_rank, ierr )

523

PEXSI_worker = (spatial_rank < npPerPole) ! Could be spin up or spin down

524

525

! PEXSI blocksize

526

pbs = norbs/npPerPole

527

528

! Careful with this. For the purposes of matrix transfers,

529

! we need the ranks of the Pole group

530

! in the "bridge" communicator/group (World)

531

532

allocate(pexsi_pole_ranks_in_world(npPerPole))

533

call MPI_Comm_Group(World_Comm, World_Group, Ierr)

534

535

call MPI_Group_translate_ranks( PEXSI_Pole_Group, npPerPole, &

536

(/ (i,i=0,npPerPole-1) /), &

537

World_Group, pexsi_pole_ranks_in_world, ierr )

538

539

! What we need is to include the actual world ranks

540

! in the distribution object

541

allocate (PEXSI_Pole_ranks_in_World_Spin(npPerPole,nspin))

542

call MPI_AllGather(pexsi_pole_ranks_in_world,npPerPole,MPI_integer,&

543

PEXSI_Pole_Ranks_in_World_Spin,npPerPole, &

544

MPI_integer,PEXSI_Spin_Comm,ierr)

545

546

! Create distributions known to all nodes

547

allocate(dist2_spin(nspin))

548

do ispin = 1, nspin

549

call newDistribution(dist2_spin(ispin), World_Comm, &

550

PEXSI_Pole_Ranks_in_World_Spin(:,ispin), &

551

TYPE_PEXSI, pbs, "px dist")

552

enddo

553

deallocate(pexsi_pole_ranks_in_world,PEXSI_Pole_Ranks_in_World_Spin)

554

call MPI_Barrier(World_Comm,ierr)

555

</pre>

556

</div>

557

</div>

558

</div>

559

560

561

<h3 id="sec-3-5">3.5 Re-distribute matrices</h3>

562

563

564

This is slightly unseemly, but it works. The <code>aux_matrix</code> derived

565

types are used to store and retrieve the matrices in either side. The

566

code is in external auxiliary modules.

567

568

569

570

571

<pre class="src src-f90">type(aux_matrix), allocatable, target :: m1_spin(:)

572

type(aux_matrix) :: m2

573

type(aux_matrix), pointer :: m1

574

integer :: nrows, nnz, nnzLocal, numColLocal

575

integer, pointer, dimension(:) :: colptrLocal=> null(), rowindLocal=>null()

576

!

577

real(dp), pointer, dimension(:) :: &

578

HnzvalLocal=>null(), SnzvalLocal=>null(), &

579

DMnzvalLocal => null()

580

!

581

integer :: ispin, pexsi_spin

582

</pre>

583

</div>

584

585

586

<pre class="src src-f90">pexsi_spin = spin_rank+1 ! {1,2}

587

! This is done serially on the Siesta side, each time

588

! filling in the structures in one PEXSI set

589

590

allocate(m1_spin(nspin))

591

do ispin = 1, nspin

592

593

m1 => m1_spin(ispin)

594

595

if (SIESTA_worker) then

596

m1%norbs = norbs

597

m1%no_l = no_l

598

m1%nnzl = sum(numH(1:no_l))

599

m1%numcols => numH

600

m1%cols => listH

601

allocate(m1%vals(2))

602

m1%vals(1)%data => S(:)

603

m1%vals(2)%data => H(:,ispin)

604

605

endif ! SIESTA_worker

606

607

call timer("redist_orbs_fwd", 1)

608

609

! Note that we cannot simply wrap this in a pexsi_spin test, as

610

! there are Siesta nodes in both spin sets.

611

! We must discriminate the PEXSI workers by the distribution info

612

dist2 => dist2_spin(ispin)

613

call redistribute_spmatrix(norbs,m1,dist1,m2,dist2,World_Comm)

614

615

call timer("redist_orbs_fwd", 2)

616

617

if (PEXSI_worker .and. (pexsi_spin == ispin) ) then

618

619

nrows = m2%norbs ! or simply 'norbs'

620

numColLocal = m2%no_l

621

nnzLocal = m2%nnzl

622

call MPI_AllReduce(nnzLocal,nnz,1,MPI_integer,MPI_sum,PEXSI_Pole_Comm,ierr)

623

624

call re_alloc(colptrLocal,1,numColLocal+1,"colptrLocal","pexsi_ldos")

625

colptrLocal(1) = 1

626

do ih = 1,numColLocal

627

colptrLocal(ih+1) = colptrLocal(ih) + m2%numcols(ih)

628

enddo

629

630

rowindLocal => m2%cols

631

SnzvalLocal => m2%vals(1)%data

632

HnzvalLocal => m2%vals(2)%data

633

634

call re_alloc(DMnzvalLocal,1,nnzLocal,"DMnzvalLocal","pexsi_ldos")

635

636

call memory_all("after setting up H+S for PEXSI (PEXSI_workers)",PEXSI_Pole_Comm)

637

638

endif ! PEXSI worker

639

enddo

640

641

! Make these available to all

642

! (Note that the values are those on process 0, which is in the spin=1 set

643

! In fact, they are only needed for calls to the interface, so the broadcast

644

! could be over PEXSI_Spatial_Comm only.

645

646

call MPI_Bcast(nrows,1,MPI_integer,0,World_Comm,ierr)

647

call MPI_Bcast(nnz,1,MPI_integer,0,World_Comm,ierr)

648

649

call memory_all("after setting up H+S for PEXSI",World_comm)

650

</pre>

651

</div>

652

</div>

653

</div>

654

655

656

<h3 id="sec-3-6">3.6 Set options</h3>

657

658

659

We use the options interface to get a template with default values,

660

and then fill in a few custom options based on fdf variables.

661

662

663

664

665

<pre class="src src-f90">type(f_ppexsi_options) :: options

666

!

667

integer :: isSIdentity

668

integer :: verbosity

669

</pre>

670

</div>

671

672

673

674

<pre class="src src-f90">! --

675

isSIdentity = 0

676

!

677

call f_ppexsi_set_default_options( options )

678

! Ordering flag:

679

! 1: Use METIS

680

! 0: Use PARMETIS/PTSCOTCH

681

options%ordering = fdf_get("PEXSI.ordering",1)

682

! Number of processors for symbolic factorization

683

! Only relevant for PARMETIS/PT_SCOTCH

684

options%npSymbFact = fdf_get("PEXSI.np-symbfact",1)

685

options%verbosity = fdf_get("PEXSI.verbosity",1)

686

</pre>

687

</div>

688

</div>

689

</div>

690

691

692

<h3 id="sec-3-7">3.7 Prepare plan</h3>

693

694

695

Each spin-set of PEXSI processors has its own plan, but we only

696

include the first-pole group of nodes…

697

698

699

700

<pre class="src src-f90">integer(c_intptr_t) :: plan

701

integer :: numProcRow, numProcCol

702

integer :: outputFileIndex

703

</pre>

704

</div>

705

706

707

708

<pre class="src src-f90"> call get_row_col(npPerPole,numProcRow,numProcCol)

709

710

! Set the outputFileIndex to be the pole index.

711

! Starting from PEXSI v0.8.0, the first processor for each pole outputs

712

! information

713

714

if( mod( mpirank, npPerPole ) .eq. 0 ) then

715

outputFileIndex = mpirank / npPerPole;

716

else

717

outputFileIndex = -1;

718

endif

719

!

720

! Note that we only use one pole's worth of processors in

721

! the global PEXSI communicator

722

!

723

plan = f_ppexsi_plan_initialize(&

724

PEXSI_Pole_Comm,&

725

numProcRow,&

726

numProcCol,&

727

outputFileIndex,&

728

info)

729

730

call check_info(info,"plan_initialize in LDOS")

731

</pre>

732

</div>

733

</div>

734

</div>

735

736

737

<h3 id="sec-3-8">3.8 Load H and S matrices</h3>

738

739

740

In this version H and S are symmetric. We associate them with the plan

741

(I really do not know very well what happens behind the

742

scenes. Presumably no copy is made.)

743

744

745

746

747

<pre class="src src-f90">call f_ppexsi_load_real_symmetric_hs_matrix(&

748

plan,&

749

options,&

750

nrows,&

751

nnz,&

752

nnzLocal,&

753

numColLocal,&

754

colptrLocal,&

755

rowindLocal,&

756

HnzvalLocal,&

757

isSIdentity,&

758

SnzvalLocal,&

759

info)

760

761

call check_info(info,"load_real_sym_hs_matrix in LDOS")

762

</pre>

763

</div>

764

</div>

765

</div>

766

767

768

<h3 id="sec-3-9">3.9 Factorization</h3>

769

770

771

This is a bit ambiguous, as we have loaded a "symmetric" matrix

772

(actually H and S), but I believe that inside (and in the plan)

773

specifically complex structures are handled and filled in.

774

775

776

777

778

<pre class="src src-f90"> call f_ppexsi_symbolic_factorize_complex_symmetric_matrix(&

779

plan, &

780

options,&

781

info)

782

call check_info(info,"factorize complex matrix in LDOS")

783

784

options%isSymbolicFactorize = 0 ! We do not need it anymore

785

</pre>

786

</div>

787

</div>

788

</div>

789

790

791

<h3 id="sec-3-10">3.10 Compute the LDOS</h3>

792

793

794

795

<pre class="src src-f90">real(dp), pointer, dimension(:) :: AnzvalLocal => null()

796

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

797

integer :: loc

798

</pre>

799

</div>

800

801

802

Note that only the first-pole team does this.

803

804

805

806

807

<pre class="src src-f90">if (PEXSI_worker) then

808

809

if(mpirank == 0) then

810

write(6,"(a,f16.5,f10.5)") &

811

'Calling PEXSI LDOS routine. Energy and broadening (eV) ', &

812

energy/eV, broadening/eV

813

write(6,"(a,i4)") &

814

'Processors working on selected inversion: ', npPerPole

815

endif

816

817

call timer("pexsi-ldos-selinv", 1)

818

819

! Build AnzvalLocal as H-zS, where z=(E,broadening), and treat

820

! it as a one-dimensional real array with 2*nnzlocal entries

821

822

call re_alloc(AnzvalLocal,1,2*nnzLocal,"AnzvalLocal","pexsi_ldos")

823

call re_alloc(AinvnzvalLocal,1,2*nnzLocal,"AinvnzvalLocal","pexsi_ldos")

824

825

loc = 1

826

do i = 1, nnzLocal

827

AnzvalLocal(loc) = Hnzvallocal(i) - energy*Snzvallocal(i)

828

AnzvalLocal(loc+1) = - broadening*Snzvallocal(i)

829

loc = loc + 2

830

enddo

831

832

call f_ppexsi_selinv_complex_symmetric_matrix(&

833

plan,&

834

options,&

835

AnzvalLocal,&

836

AinvnzvalLocal,&

837

info)

838

839

call check_info(info,"selinv complex matrix in LDOS")

840

841

! Get DMnzvalLocal as 1/pi * Imag(Ainv...)

842

843

loc = 1

844

do i = 1, nnzLocal

845

DMnzvalLocal(i) = (1.0_dp/pi) * AinvnzvalLocal(loc+1)

846

loc = loc + 2

847

enddo

848

call de_alloc(AnzvalLocal,"AnzvalLocal","pexsi_ldos")

849

call de_alloc(AinvnzvalLocal,"AinvnzvalLocal","pexsi_ldos")

850

851

call timer("pexsi-ldos-selinv", 2)

852

!

853

endif ! PEXSI_worker

854

</pre>

855

</div>

856

</div>

857

</div>

858

859

860

<h3 id="sec-3-11">3.11 Copy information to Siesta side</h3>

861

862

863

864

<pre class="src src-f90">do ispin = 1, nspin

865

866

m1 => m1_spin(ispin)

867

868

if (PEXSI_worker .and. (pexsi_spin == ispin)) then

869

! Prepare m2 to transfer

870

871

call de_alloc(colPtrLocal,"colPtrLocal","pexsi_ldos")

872

873

call de_alloc(m2%vals(1)%data,"m2%vals(1)%data","pexsi_ldos")

874

call de_alloc(m2%vals(2)%data,"m2%vals(2)%data","pexsi_ldos")

875

876

deallocate(m2%vals)

877

allocate(m2%vals(1))

878

m2%vals(1)%data => DMnzvalLocal(1:nnzLocal)

879

880

endif

881

882

! Prepare m1 to receive the results

883

if (SIESTA_worker) then

884

nullify(m1%vals(1)%data) ! formerly pointing to S

885

nullify(m1%vals(2)%data) ! formerly pointing to H

886

deallocate(m1%vals)

887

nullify(m1%numcols) ! formerly pointing to numH

888

nullify(m1%cols) ! formerly pointing to listH

889

endif

890

891

call timer("redist_orbs_bck", 1)

892

dist2 => dist2_spin(ispin)

893

call redistribute_spmatrix(norbs,m2,dist2,m1,dist1,World_Comm)

894

call timer("redist_orbs_bck", 2)

895

896

if (PEXSI_worker .and. (pexsi_spin == ispin)) then

897

call de_alloc(DMnzvalLocal, "DMnzvalLocal", "pexsi_ldos")

898

899

nullify(m2%vals(1)%data) ! formerly pointing to DM

900

deallocate(m2%vals)

901

! allocated in the direct transfer

902

call de_alloc(m2%numcols,"m2%numcols","pexsi_ldos")

903

call de_alloc(m2%cols, "m2%cols", "pexsi_ldos")

904

endif

905

906

if (SIESTA_worker) then

907

908

LocalDOSDM(:,ispin) = m1%vals(1)%data(:)

909

! Check no_l

910

if (no_l /= m1%no_l) then

911

call die("Mismatch in no_l")

912

endif

913

! Check listH

914

if (any(listH(:) /= m1%cols(:))) then

915

call die("Mismatch in listH")

916

endif

917

918

call de_alloc(m1%vals(1)%data,"m1%vals(1)%data","pexsi_ldos")

919

deallocate(m1%vals)

920

! allocated in the direct transfer

921

call de_alloc(m1%numcols,"m1%numcols","pexsi_ldos")

922

call de_alloc(m1%cols, "m1%cols", "pexsi_ldos")

923

924

endif

925

enddo

926

call timer("pexsi-ldos", 2)

927

</pre>

928

</div>

929

</div>

930

</div>

931

932

933

<h3 id="sec-3-12">3.12 Clean up</h3>

934

935

936

937

<pre class="src src-f90">call delete(dist1)

938

do ispin = 1, nspin

939

call delete(dist2_spin(ispin))

940

enddo

941

deallocate(dist2_spin)

942

deallocate(m1_spin)

943

944

call MPI_Comm_Free(PEXSI_Spatial_Comm, ierr)

945

call MPI_Comm_Free(PEXSI_Spin_Comm, ierr)

946

call MPI_Comm_Free(World_Comm, ierr)

947

948

! This communicator was created from a subgroup.

949

! As such, it is MPI_COMM_NULL for those processes

950

! not in the subgroup (non PEXSI_workers). Only

951

! defined communicators can be freed

952

! Also, we finalize the plan here

953

if (PEXSI_worker) then

954

call f_ppexsi_plan_finalize( plan, info )

955

call MPI_Comm_Free(PEXSI_Pole_Comm, ierr)

956

endif

957

958

call MPI_Group_Free(PEXSI_Spatial_Group, ierr)

959

call MPI_Group_Free(PEXSI_Pole_Group, ierr)

960

call MPI_Group_Free(World_Group, ierr)

961

</pre>

962

</div>

963

</div>

964

</div>

965

966

967

<h3 id="sec-3-13">3.13 Support routines</h3>

968

969

970

A couple of routines

971

972

973

974

975

976

<<check-info>>

977

</pre>

978

</div>

979

</div>

980

981

982

<h4 id="sec-3-13-1">3.13.1 Row and column partition of npPerPole</h4>

983

984

985

986

<pre class="src src-f90">subroutine get_row_col(np,nrow,ncol)

987

integer, intent(in) :: np

988

integer, intent(out) :: nrow, ncol

989

!

990

! Finds the factors nrow and ncol such that nrow*ncol=np,

991

! are as similar as possible, and nrow>=ncol.

992

! For prime np, ncol=1, nrow=np.

993

994

ncol = floor(sqrt(dble(np)))

995

do

996

nrow = np/ncol

997

if (nrow*ncol == np) exit

998

ncol = ncol - 1

999

enddo

1000

end subroutine get_row_col

1001

</pre>

1002

</div>

1003

</div>

1004

</div>

1005

1006

<h4 id="sec-3-13-2">3.13.2 Error dispatcher</h4>

1007

1008

1009

1010

<pre class="src src-f90">subroutine check_info(info,str)

1011

integer, intent(in) :: info

1012

character(len=*), intent(in) :: str

1013

1014

if(mpirank == 0) then

1015

if (info /= 0) then

1016

write(6,*) trim(str) // " info : ", info

1017

call die("Error exit from " // trim(str) // " routine")

1018

endif

1019

call pxfflush(6)

1020

endif

1021

end subroutine check_info

1022

</pre>

1023

</div>

1024

</div>

1025

</div>

1026

</div>

1027

</div>

1028

</div>

1029

1030

Author: Alberto Garcia

1031

Created: 2015-11-05 Thu 12:35

1032

<a href="http://www.gnu.org/software/emacs/">Emacs</a> 24.5.1 (<a href="http://orgmode.org">Org</a> mode 8.2.10)

1033

<a href="http://validator.w3.org/check?uri=referer">Validate</a>

1034

</div>

1035

</body>

1036

</html>