1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
|
!
! Copyright (C) 1996-2016 The SIESTA group
! This file is distributed under the terms of the
! GNU General Public License: see COPYING in the top directory
! or http://www.gnu.org/copyleft/gpl.txt.
! See Docs/Contributors.txt for a list of contributors.
!
subroutine mulliken(iopt,natoms,nbasistot,maxnh,numh,
. listhptr,listh,s,dm,isa,lasto,iaorb,
. iphorb)
C ********************************************************************
C Subroutine to perform Mulliken population analysis.
C (Overlap and total populations, both for orbitals and for atoms)
C The density matrix (d.m.) and overlap matrix are passed in sparse form
C (both with the same sparse structure)
C There is no output. The populations are printed to the output.
C
C Written by P.Ordejon, October'96
C Non-collinear spin added by J.M.Soler, May 1998.
C Symmetry label for each orbital included by DSP, Oct. 1998.
C Label with the principal quantum number introduced by DSP, July 1999.
C ************************** INPUT ************************************
C integer iopt : Work option: 1 = atomic and orbital charges
C 2 = 1 + atomic overlap pop.
C 3 = 2 + orbital overlap pop.
C integer natoms : Number of atoms
C integer nbasistot : Number of basis orbitals over all nodes
C integer maxnh : First dimension of d.m. and overlap, and its
C maximum number of non-zero elements
C integer numh(nbasis) : First Control vector of d.m. and overlap
C integer listhptr(nbasis) : Second Control vector of d.m. and overlap
C integer listh(maxnh) : Third Control vector of d.m. and overlap
C real*8 s(maxnh) : Overlap matrix in sparse form
C real*8 dm(maxnh,spin%DM) : Density matrix in sparse form
C integer isa(natoms) : Species index of each atom
C integer lasto(0:maxa) : Index of last orbital of each atom
C (lasto(0) = 0)
C integer iaorb(nbasis) : Atomic index of each orbital 19684
C integer iphorb(nbasis) : Orbital index of each orbital in its atom
C ************************* OUTPUT *************************************
C No output. The results are printed to standard output
C **********************************************************************
C
C Modules
C
use precision, only : dp
use parallel, only : IOnode, Node, Nodes
use parallelsubs, only : GlobalToLocalOrb, GetNodeOrbs
use parallelsubs, only : LocalToGlobalOrb
use atmfuncs, only : symfio, cnfigfio, labelfis, nofis
use alloc, only : re_alloc, de_alloc
use siesta_cml
use m_spin, only : spin
#ifdef MPI
use mpi_siesta
#endif
implicit none
integer, intent(in) ::
. iopt,natoms,nbasistot,maxnh
integer
. numh(*),lasto(0:natoms),listh(maxnh),listhptr(*),
. iphorb(*), isa(natoms), iaorb(*)
real(dp)
. dm(maxnh,spin%DM),s(maxnh)
C Internal parameters ..................................................
C Number of columns in printout. Must be smaller than 20
integer ncol, nbasis
parameter (ncol = 8)
#ifdef MPI
integer :: MPIerror, mpistatus(MPI_STATUS_SIZE)
real(dp) :: pb(ncol)
real(dp), pointer :: qb(:)=>null()
character(len=128) :: mpibuff ! 128 is sufficiently long
#endif
integer i,ia,ib,ii,imax,in,io,ior,ip,is,j,ja,jja,jo,jor,
. ind, nao, nblock, ns, ispec, irow, nrow,
. config(ncol), itot
real(dp)
. p(ncol),qa,qas(4),qts(4),qtot,qtot_temp,stot,svec(3)
real(dp) :: dm_aux(4)
real(dp), pointer :: qo(:)=>null(), qos(:,:)=>null()
character sym_label(ncol)*7, atm_label*20
character(len=8), pointer :: orb_label(:)=>null()
character(len=13) :: atomtitle
C ......................
#ifdef DEBUG
call write_debug( ' PRE mulliken' )
#endif
C Get Node numberReduce
#ifdef MPI
C Find number of locally stored orbitals and allocated related arrays
call GetNodeOrbs(nbasistot,Node,Nodes,nbasis)
#else
nbasis = nbasistot
#endif
if (iopt.eq.0) then
C iopt = 0 implies no analysis
#ifdef DEBUG
call write_debug( ' POS mulliken iopt==0' )
#endif
return
elseif (iopt.lt.0 .or. iopt.gt.3) then
if (ionode) then
write(6,"(a)") 'mulliken: ERROR: Wrong iopt'
endif
#ifdef DEBUG
call write_debug( ' POS mulliken iopt<0 or iopt>3' )
#endif
return
endif
ns=0
do i = 1,natoms
ns=max(ns,isa(i))
enddo
C Compute and print Overlap Populations for Orbitals ....................
if (iopt .eq. 3) then
if (ionode) then
write(6,*)
write(6,"(a)")'mulliken: Overlap Populations between Orbitals'
endif
nblock = nbasistot / ncol
ip=1
if (nblock*ncol .eq. nbasistot) ip=0
do ib = 1,nblock+ip
imax = ncol
if (ib .eq. nblock+1) imax = nbasistot - nblock * ncol
do ii=1,imax
sym_label(ii)=symfio(isa(iaorb((ib-1)*ncol+ii)),
. iphorb((ib-1)*ncol+ii))
config(ii)=cnfigfio(isa(iaorb((ib-1)*ncol+ii)),
. iphorb((ib-1)*ncol+ii))
enddo
if (ionode) then
write(6,*)
write(6,'(14x,20(2x,i4,3x))')((ib-1)*ncol+ii,ii=1,imax)
write(6,'(17x,20(1x,i1,a7))')
. (config(ii),sym_label(ii),ii=1,imax)
endif
! Following algorithm is designed to ensure that terms relating
! to orbitals are printed out in the correct sequence. It's even
! less efficient than the algorithm that used to be here - on the
! other hand, it does actually work properly.
! (Briefly: iterate over orbitals -
! if it's on the current node, and we are node 0 then print immediately;
! else if it's not, print the output to a buffer & send it to node 0.
! if it's not on the current node, then if we're node 0, wait for the
! message from whichever node it is on, and print the buffer.
! Otherwise do nothing.
do itot = 1,nbasistot
call GlobalToLocalOrb(itot,Node,Nodes,i)
if (i.gt.0) then
sym_label(1)=symfio(isa(iaorb(itot)),iphorb(itot))
config(1)=cnfigfio(isa(iaorb(itot)),iphorb(itot))
do ii = 1,imax
p(ii) = 0._dp
enddo
do in = 1,numh(i)
ind = listhptr(i)+in
j = listh(ind)
ii = j - (ib - 1) * ncol
if (ii .ge. 1 .and. ii .le. imax) then
p(ii) = 0._dp
do is = 1,spin%Spinor
p(ii) = p(ii) + dm(ind,is) * s(ind)
enddo
endif
enddo
if (ionode) then
write(6,15) itot,config(1), sym_label(1),
. (p(ii),ii=1,imax)
#ifdef MPI
else
write(mpibuff,15) itot,config(1), sym_label(1),
. (p(ii),ii=1,imax)
call mpi_send(mpibuff, 128, MPI_Character, 0, itot,
. MPI_Comm_world, mpierror)
endif
else
if (ionode) then
call mpi_recv(mpibuff, 128, MPI_Character,
. MPI_ANY_SOURCE, itot,
. MPI_Comm_world, mpistatus, mpierror)
write(6,'(a)') trim(mpibuff)
#endif
endif
endif
enddo
enddo
endif
C ...................
C Compute and print Overlap Populations for Atoms ....................
if (iopt .ge. 2) then
if (ionode) then
write(6,*)
write(6,"(a)")'mulliken: Overlap Populations between Atoms'
endif
nblock = natoms / ncol
ip=1
if (nblock*ncol .eq. natoms) ip=0
do ib = 1,nblock+ip
imax = ncol
if (ib .eq. nblock+1) imax = natoms - nblock * ncol
if (ionode) then
write(6,*)
write(6,10) ((ib-1)*ncol+ii,ii=1,imax)
endif
do i = 1,natoms
do ii = 1,imax
p(ii) = 0.0
enddo
do ior = lasto(i-1)+1,lasto(i)
call GlobalToLocalOrb(ior,Node,Nodes,itot)
if (itot.gt.0) then
do in = 1,numh(itot)
ind = listhptr(itot)+in
jor = listh(ind)
do jja = 1,natoms
if (lasto(jja) .ge. jor) then
ja = jja
goto 100
endif
enddo
goto 110
100 ii = ja - (ib - 1) * ncol
if (ii .ge. 1 .and. ii .le. imax) then
do is = 1,spin%Spinor
p(ii) = p(ii) + dm(ind,is) * s(ind)
enddo
endif
110 continue
enddo
endif
enddo
#ifdef MPI
C Global sum of values stored in p
pb(1:imax)=p(1:imax)
call MPI_Reduce(pb,p,imax,MPI_double_precision,MPI_sum,
. 0,MPI_Comm_World,MPIerror)
#endif
if (ionode) then
write(6,11) i,(p(ii),ii=1,imax)
endif
enddo
enddo
endif
C ....................
C Compute and print Mulliken Orbital and Atomic Populations ..........
if (iopt .ge. 1) then
if (ionode) then
write(6,*)
write(6,"(a)")'mulliken: Atomic and Orbital Populations:'
endif
if (spin%DM .le. 2) then
C We only write mulliken to cml file if we are unpolarized or collinear
C and (currently) don't bother in the non-collinear case.
if (cml_p) then
call cmlStartPropertyList(xf=mainXML,
. dictRef='siesta:mulliken',
. title='Mulliken Population Analysis')
endif
do is = 1,spin%DM
if (spin%DM .eq. 2) then
if(is .eq. 1) then
if (ionode) write(6,'(/,a)') 'mulliken: Spin UP '
if (cml_p) call cmlStartPropertyList(xf=mainXML,
. title='SpinUp')
elseif(is .eq. 2) then
if (ionode) write(6,'(/,a)') 'mulliken: Spin DOWN '
if (cml_p) call cmlStartPropertyList(xf=mainXML,
. title='SpinDown')
endif
else
if (cml_p) call cmlStartPropertyList(xf=mainXML,
. title='Unpolarized')
endif
qtot = 0._dp
do ispec = 1,ns
atm_label = labelfis(ispec)
if (ionode) then
write(6,'(/2a)')'Species: ', atm_label
write(6,'(a4,a7,a6)') 'Atom', 'Qatom', 'Qorb'
endif
C DSP, Writing symmetries for each orbital.
C DSP, Orbitals with a 'P' belong to the polarization shell
nao = nofis(ispec)
call re_alloc( orb_label, 1, nao, 'orb_label', 'mulliken' )
do io=1,nao
write(orb_label(io),'(i1,a7)')
. cnfigfio(ispec,io), symfio(ispec,io)
enddo
nrow=nao/ncol
io=1
do irow=1,nrow
if (ionode) then
write(6,'(15x,8(a8))') orb_label(io:io+ncol-1)
io = io+ncol
endif
enddo
if (ionode) then
write(6,'(15x,8(a8))') orb_label(io:nao)
endif
if (cml_p) then
call cmlStartPropertyList(xf=mainXML,
. title=trim(atm_label))
call cmlAddProperty(xf=mainXML,
. title='Orbital list', value=orb_label)
endif
call de_alloc( orb_label, 'orb_label', 'mulliken' )
do ia = 1,natoms
if (isa(ia).eq.ispec) then
C Compute charge in each orbital of atom ia
nao = lasto(ia) - lasto(ia-1)
call re_alloc( qo, 1, nao, 'qo', 'mulliken' )
qa = 0.0_dp
do io = lasto(ia-1)+1,lasto(ia)
nao = io - lasto(ia-1)
qo(nao) = 0.0_dp
call GlobalToLocalOrb(io,Node,Nodes,itot)
if (itot.gt.0) then
do in = 1,numh(itot)
ind = listhptr(itot) + in
qo(nao) = qo(nao) + dm(ind,is) * s(ind)
enddo
qa = qa + qo(nao)
endif
enddo
qtot = qtot + qa
#ifdef MPI
C Global sum of values stored in qo
call re_alloc( qb, 1, nao, 'qb', 'mulliken' )
qb = qo
call MPI_Reduce(qb,qo,nao,MPI_double_precision,
. MPI_sum,0,MPI_Comm_World,MPIerror)
qb(1)=qa
call MPI_Reduce(qb(1),qa,1,MPI_double_precision,
. MPI_sum,0,MPI_Comm_World,MPIerror)
call de_alloc( qb, 'qb', 'mulliken' )
#endif
if (ionode) then
write(6,'(i4,f7.3,8f8.3,(/11x,8f8.3))')
. ia, qa, qo
endif
if (cml_p) then
write(atomtitle,'(a,i8)') 'atom ',ia
call cmlStartPropertyList(xf=mainXML,
. title=trim(atomtitle))
call cmlAddProperty(xf=mainXML,
. title="Atomic Charge", value=qa,
. units="siestaUnits:e")
call cmlAddProperty(xf=mainXML,
. title="Orbital Charges", value=qo,
. units="siestaUnits:e")
call cmlEndPropertyList(xf=mainXML) !atomtitle
endif
call de_alloc( qo, 'qo', 'mulliken' )
endif
enddo
if (cml_p) call cmlEndPropertyList(xf=mainXML) !atm_label
enddo
#ifdef MPI
C Global sum of total charge
qtot_temp = qtot
call MPI_Reduce(qtot_temp,qtot,1,MPI_double_precision,
$ MPI_sum, 0,MPI_Comm_World,MPIerror)
#endif
if (ionode) then
write(6,"(/a,f12.3)") 'mulliken: Qtot = ', qtot
endif
if (cml_p) then
call cmlAddProperty(xf=mainXML,
. title="Total Charge", value=qtot,
. units="siestaUnits:e")
call cmlEndPropertyList(xf=mainXML) ! spins
endif
enddo
if (cml_p) call cmlEndPropertyList(xf=mainXML) ! mulliken
elseif (spin%DM .ge. 4) then
call re_alloc( qos, 1, spin%Grid, 1, nbasistot,
& 'qos', 'mulliken' )
do is = 1,spin%Grid
qts(is) = 0.0_dp
do io = 1,nbasistot
qos(is,io) = 0.0_dp
enddo
enddo
! For SOC-onsite, dm(:,4) has the wrong sign, but this
! routine is not called ('moments' is).
do io = 1,nbasis
call LocalToGlobalOrb(io,Node, Nodes, itot)
do in = 1,numh(io)
ind = listhptr(io)+in
jo = listh(ind)
dm_aux(1:4) = dm(ind,1:4)
if (spin%SO) then
dm_aux(3) = 0.5_dp*(dm_aux(3)+dm(ind,7))
dm_aux(4) = 0.5_dp*(dm_aux(4)+dm(ind,8))
endif
do is = 1,spin%Grid
qos(is,itot) = qos(is,itot) + dm_aux(is)*s(ind)/2
if ( jo <= nbasistot )
. qos(is,jo) = qos(is,jo) + dm_aux(is)*s(ind)/2
enddo
enddo
enddo
#ifdef MPI
call re_alloc( qb, 1, spin%Grid*nbasistot,
& 'qb', 'mulliken' )
call MPI_Reduce(qos(1,1), qb(1), spin%Grid*nbasistot,
& Mpi_Double_Precision, MPI_Sum,
& 0, MPI_Comm_world, MPIerror)
itot = 0
do io = 1 , nbasistot
do is = 1 , spin%Grid
itot = itot + 1
qos(is, io) = qb(itot)
end do
end do
call de_alloc(qb, 'qb', 'mulliken')
#endif
! The above reduction means that the
! following mulliken calculation actually is
! serial.
! We could have performed other analysis, but
! that would prove additional "small" communication.
! This should aptly provide speed and flexibility.
! Also, this entire loop may be encapsulated
! In (if ionode)... However, this will also
! streamline the processors.
do ispec=1,ns
atm_label=labelfis(ispec)
if (ionode) then
write(6,'(/2a)')'Species: ', atm_label
write(6,'(/,a4,1x,a8,4x,2a10,3x,a8,/,a)')
. 'Atom', 'Orb', 'Charge', 'Spin', 'Svec',
. repeat('-', 64)
endif
do ia = 1,natoms
! Check for same species
if ( isa(ia) /= ispec ) cycle
! Initialize atomic charges (per spin)
qas(1:spin%Grid) = 0._dp
do io = lasto(ia-1)+1, lasto(ia)
C DSP, Writing symmetries for each orbital.
C DSP, Orbitals with a 'P' belong to the polarization shell
sym_label(1) = symfio(ispec,iphorb(io))
config(1) = cnfigfio(ispec,iphorb(io))
! Reduce total charge per-atom
do is = 1 , spin%Grid
qas(is) = qas(is) + qos(is,io)
enddo
call spnvec( spin%Grid, qos(1:spin%Grid,io),
$ qtot, stot, svec)
if ( IONode ) then
write(6,'(i5,i3,tr1,i1,a7,2f10.5,3x,3f8.3)')
. ia,io-lasto(ia-1),config(1),sym_label(1),
. qtot, stot, svec
end if
end do ! atomic orbitals
call spnvec( spin%Grid, qas, qtot, stot, svec )
if ( IONode ) then
write(6,'(i5,5x,a5,2x,2f10.5,3x,3f8.3,/)')
. ia, 'Total', qtot, stot, svec
end if
! Reduce total charge in system
do is = 1 , spin%Grid
qts(is) = qts(is) + qas(is)
end do
end do ! loop atoms
call spnvec( spin%Grid, qts, qtot, stot, svec )
if ( IONode ) then
write(6,'(a,/,a5,12x,2f10.5,3x,3f8.3,/)')
. repeat('-',64), 'Total', qtot, stot, svec
end if
end do ! species
call de_alloc( qos, 'qos', 'mulliken' )
end if ! spin%DM .ge. 4
end if
C ...................
#ifdef DEBUG
call write_debug( ' POS mulliken' )
#endif
10 format(12x,20(2x,i4,2x))
11 format(i12,20(1x,f7.3))
15 format(i4,1x,i1,a7,20(1x,f8.3))
return
end subroutine mulliken
subroutine spnvec( ns, qs, qt, DPOL, sv )
c ********************************************************************
c Finds the spin vector components from the spin density matrix
c Written by J.M.Soler, May 1998.
c ******* Input ******************************************************
c integer ns : Number of components in spin density matrix
c real*8 qs(ns) : Spin density matrix elements with the convention
c is=1 => Q11; is=2 => Q22; is=3 => Real(Q12);
c is=4 => Imag(Q12)
c ******* Output *****************************************************
c real*8 qt : Total charge
c real*8 DPOL : Total spin
c real*8 sv(3) : Spin vector
c ********************************************************************
use precision
implicit none
integer, intent(in) :: ns
real(dp), intent(in) :: qs(ns)
real(dp), intent(out) :: qt, DPOL, sv(3)
! Local variables
real(dp) :: dxy, dz, cosph, costh, sinph, sinth
real(dp), parameter :: DENMIN = 1.e-12_dp
if (ns .eq. 1) then
qt = qs(1)
DPOL = 0.0_dp
sv(1) = 0.0_dp
sv(2) = 0.0_dp
sv(3) = 0.0_dp
elseif (ns .eq. 2) then
qt = qs(1) + qs(2)
DPOL = qs(1) - qs(2)
sv(1) = 0.0_dp
sv(2) = 0.0_dp
sv(3) = DPOL
elseif (ns .eq. 4) then
! TODO, correct for possible division by zero (instead of DENMIN)
qt = qs(1) + qs(2)
dz = qs(1) - qs(2)
dxy = 2._dp * sqrt( qs(3) ** 2 + qs(4) ** 2 )
DPOL = sqrt( dz ** 2 + dxy ** 2 )
DPOL = max( DPOL , DENMIN )
costh = dz / DPOL
sinth = sqrt( 1.0_dp - costh**2 )
dxy = max( dxy , DENMIN ) * .5_dp
cosph = qs(3) / dxy
sinph = qs(4) / dxy
sv(1) = DPOL * sinth * cosph
sv(2) = DPOL * sinth * sinph
sv(3) = DPOL * costh
else
write(6,*) 'spnvec: ERROR: invalid argument ns =', ns
return
endif
end subroutine spnvec
|