← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Util/COOP/mprop.f90

  • Committer: Alberto Garcia
  • Date: 2019-09-02 14:09:43 UTC
  • mfrom: (427.6.323 trunk)
  • Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Util/COOP/mprop.f90
17
17
 
18
18
  implicit none
19
19
 
20
 
  logical :: gamma_wfsx, got_qcos
 
20
  logical :: gamma_wfsx, got_qcos, non_coll
21
21
  integer :: ii1, ii2, ind, ind_red, no1, no2, n_int, nnz
22
 
  integer :: imin, imax
23
 
  real(dp) :: factor
 
22
  integer :: imin, imax, wfs_spin_flag, nspin_blocks
 
23
  
 
24
  complex(DP), dimension(2)    :: spinor_1, spinor_2, H_c2
 
25
  complex(DP), dimension(2,2)  :: H
 
26
  complex(DP) :: c_c1_c2, c_c1_H_c2
 
27
  real(dp) :: factor_S, factor_H, dos_value
24
28
 
25
29
  ! We use a smearing function of the form f(x) = exp(-(x/smear)**2) / (smear*sqrt(pi))
26
30
  ! A weight tolerance of 1.0e-4 corresponds to going about 3*smear on either
136
140
  read(wfs_u) nkp, gamma_wfsx
137
141
  allocate (wk(nkp), pk(3,nkp))
138
142
 
139
 
  read(wfs_u) nsp
 
143
  read(wfs_u) wfs_spin_flag   !  1, 2, or 4 
 
144
  non_coll = (wfs_spin_flag == 4)
140
145
  read(wfs_u) nao
141
146
  read(wfs_u)        !! Symbols, etc
142
 
  if (debug) print *, "WFSX read: nkp, nsp, nnao: ", nkp, nsp, nao
143
 
 
144
 
  allocate (ados(npts_energy,nsp), ww(npts_energy))
 
147
  if (debug) print *, "WFSX read: nkp, nspin_flag, nnao: ", nkp, wfs_spin_flag, nao
 
148
 
 
149
  if (non_coll) then
 
150
     nspin_blocks = 1
 
151
  else
 
152
     nspin_blocks = wfs_spin_flag
 
153
  endif
 
154
  print *, "NSPIN BLOCKS: ", nspin_blocks
 
155
 
 
156
  allocate (ados(npts_energy,nspin_blocks), ww(npts_energy))
145
157
 
146
158
  nwfmx = -huge(1)
147
159
  nwfmin = huge(1)
150
162
  min_eigval_in_band_set = huge(1.0_dp)
151
163
  max_eigval_in_band_set = -huge(1.0_dp)
152
164
 
 
165
 
153
166
  do ik=1,nkp
154
 
     do is=1,nsp
 
167
     do is=1,nspin_blocks
155
168
 
156
169
        read(wfs_u) idummy, pk(1:3,ik), wk(ik)
157
170
        if (idummy /= ik) stop "ik index mismatch in WFS file"
252
265
 
253
266
 
254
267
  e_step = (high_e-low_e)/(npts_energy-1)
255
 
  ados(:,1:nsp) = 0.0_dp
 
268
  ados(:,1:nspin_blocks) = 0.0_dp
256
269
 
257
270
  ! skip four records
258
271
 
264
277
  read(wfs_u) 
265
278
 
266
279
  do ik=1,nkp
267
 
     do is=1,nsp
 
280
     do is=1,nspin_blocks
268
281
        read(wfs_u)
269
282
        read(wfs_u)
270
283
        read(wfs_u)  number_of_wfns
284
297
     enddo
285
298
  enddo
286
299
 
287
 
  call io_assign(idos)
288
 
  open(idos,file=trim(sflnm)//".alldos",form="formatted", &
289
 
       status="unknown",action="write",position="rewind")
290
 
  write(idos,*) "#  Energy   LARGE-SCALE DOS"
291
 
 
292
 
  do i = 1, npts_energy
293
 
     energy = low_e + e_step*(i-1)
294
 
     write(idos,*) energy, (ados(i,is), is=1,nsp)
295
 
  enddo
296
 
 
297
 
  call io_close(idos)
 
300
  ! Write "LARGE-SCALE DOS"  
 
301
  call write_curve_to_file(trim(sflnm)//".alldos",ados)
298
302
 
299
303
  ! Read HSX file
300
304
  ! Will pick up atoms, zval, and thus the nominal number of electrons,
312
316
 
313
317
  !
314
318
  !       Compute integrated total DOS
315
 
  !       Here we double the DOS for the case of nspin=1,
316
 
  !       to get the correct number of states.
317
319
 
318
320
  allocate(intdos(npts_energy), intebs(npts_energy))
319
321
  call io_assign(intdos_u)
324
326
  write(intdos_u,*) low_e, intdos(1), intebs(1)
325
327
  do i = 2, npts_energy
326
328
     energy = low_e + e_step*(i-1)
327
 
     intdos(i) = intdos(i-1) + sum(ados(i,:)) * e_step * 2.0_dp /nsp
328
 
     intebs(i) = intebs(i-1) + energy*sum(ados(i,:)) * e_step * 2.0_dp /nsp
 
329
 
 
330
     ! For the spinless case, we need a factor of two for proper normalization here
 
331
     dos_value = sum(ados(i,:))
 
332
     if (wfs_spin_flag == 1) dos_value = 2*dos_value
 
333
 
 
334
     intdos(i) = intdos(i-1) + dos_value * e_step 
 
335
     intebs(i) = intebs(i-1) + energy*dos_value * e_step 
329
336
     write(intdos_u,*) energy, intdos(i), intebs(i)
330
337
  enddo
331
338
  call io_close(intdos_u)
500
507
  !=====================
501
508
 
502
509
 if (coop) then
503
 
    allocate (coop_vals(npts_energy,nsp,ncb))
504
 
    allocate (cohp_vals(npts_energy,nsp,ncb))
505
 
    coop_vals(:,:nsp,:ncb)=0.d0
506
 
    cohp_vals(:,:nsp,:ncb)=0.d0
 
510
    allocate (coop_vals(npts_energy,nspin_blocks,ncb))
 
511
    allocate (cohp_vals(npts_energy,nspin_blocks,ncb))
 
512
    coop_vals(:,:,:)=0.d0
 
513
    cohp_vals(:,:,:)=0.d0
507
514
 endif
508
515
 if (dos) then
509
 
    allocate (pdos_vals(npts_energy,nsp,ncb))
510
 
    pdos_vals(:,:nsp,:ncb) = 0.0_dp
 
516
    allocate (pdos_vals(npts_energy,nspin_blocks,ncb))
 
517
    pdos_vals(:,:,:) = 0.0_dp
511
518
 endif
512
519
 
513
 
  ados(:,:nsp) = 0.0_dp
 
520
  ados(:,:) = 0.0_dp
514
521
 
515
522
  !================================================================
516
523
 
517
524
  ! * Curves
518
525
 
519
 
     if (gamma_wfsx) then
520
 
        allocate(wf(1,1:no_u))
 
526
     ! The first dimension is the number of real numbers per orbital
 
527
     ! 1 for real wfs, 2 for complex, and four for the two spinor components
 
528
  
 
529
     if (non_coll) then
 
530
        allocate(wf_single(4,1:no_u))
 
531
        allocate(wf(4,1:no_u))
521
532
     else
522
 
        allocate(wf(2,1:no_u))
 
533
        if (gamma_wfsx) then
 
534
           allocate(wf_single(1,1:no_u))
 
535
           allocate(wf(1,1:no_u))
 
536
        else
 
537
           allocate(wf_single(2,1:no_u))
 
538
           allocate(wf(2,1:no_u))
 
539
        endif
523
540
     endif
524
 
 
525
541
     allocate (mask2(1:no_u))
526
542
 
527
543
     do ic=1,ncb
613
629
        if (debug) print *, "Number of k-points, spins: ", nkp, nsp
614
630
        do ik=1,nkp
615
631
           if (debug) print *, "k-point: ", ik
616
 
           do is=1,nsp
 
632
           do is=1,nspin_blocks
617
633
              read(wfs_u)
618
634
              read(wfs_u)
619
635
              read(wfs_u)  number_of_wfns
635
651
                    CYCLE
636
652
                 endif
637
653
 
638
 
                 read(wfs_u) (wf(:,io), io=1,nao)
639
 
 
 
654
                 read(wfs_u) (wf_single(:,io), io=1,no_u)
 
655
                 ! Use a double precision form in what follows
 
656
                 wf(:,:) = real(wf_single(:,:), kind=dp)
640
657
 
641
658
                 ! This block will be repeated for every curve,
642
659
                 ! but we will divide by the number of curves before writing out
670
687
                         ind_red = ptr(i1)+i2
671
688
                         io2 = list_io2(ind_red)
672
689
                         ind = list_ind(ind_red)
673
 
                             
674
 
                                ! (qcos, qsin) = C_1*conjg(C_2)
675
 
                                !AG: Corrected:  (qcos, qsin) = conjg(C_1)*(C_2)
676
 
                                ! We might want to avoid recomputing this
677
 
 
678
 
                                if (gamma_wfsx) then
679
 
                                   qcos = wf(1,io1)*wf(1,io2) 
680
 
                                   qsin = 0.0_dp
 
690
 
 
691
                                ! (qcos, qsin) = conjg(C_1) * [S or H] * C_2
 
692
                                ! We might want to avoid recomputing the "S" pure wf parts
 
693
 
 
694
                                if (non_coll) then
 
695
 
 
696
                                   ! Use 'dp' to keep the wfs in double precision
 
697
                                   ! (Recall that wf() is now dp; converted right after reading)
 
698
                                   spinor_1 = [ cmplx(wf(1,io1),wf(2,io1), dp), &
 
699
                                                      cmplx(wf(3,io1),wf(4,io1), dp) ]
 
700
                                   spinor_2 = [ cmplx(wf(1,io2),wf(2,io2), dp), &
 
701
                                                cmplx(wf(3,io2),wf(4,io2), dp) ]
 
702
 
 
703
                                   ! We take the signs for 1,2 and 2,1
 
704
                                   ! from the construction of Ebs_Haux in 'compute_energies'
 
705
                                   if (h_spin_dim == 8) then
 
706
                                      H(1,1) = cmplx(Hamilt(ind,1), Hamilt(ind,5), dp)
 
707
                                      H(1,2) = cmplx(Hamilt(ind,3), -Hamilt(ind,4), dp)
 
708
                                      H(2,1) = cmplx(Hamilt(ind,7), Hamilt(ind,8), dp)
 
709
                                      H(2,2) = cmplx(Hamilt(ind,2), Hamilt(ind,6), dp)
 
710
                                   else   ! nsp=4; just non-collinear; no SOC
 
711
                                      H(1,1) = cmplx(Hamilt(ind,1), 0.0_dp, dp)
 
712
                                      H(1,2) = cmplx(Hamilt(ind,3), -Hamilt(ind,4), dp)
 
713
                                      H(2,1) = cmplx(Hamilt(ind,3), Hamilt(ind,4), dp)
 
714
                                      H(2,2) = cmplx(Hamilt(ind,2), 0.0_dp, dp)
 
715
                                   endif
 
716
 
 
717
                                   ! For DOS/COOP, take as weight the "complete spinor" product
 
718
                                   ! The dot_product is directly sum(conjg(a1)*a2) for complex arrays
 
719
                                   ! The kind returned is "the highest" of the arguments. In this case,
 
720
                                   ! it will be 'dp' if the 'spinor_X' variables are 'dp'.
 
721
                                   c_c1_c2 = dot_product(spinor_1,spinor_2)
 
722
                                   qcos= real(c_c1_c2, dp) 
 
723
                                   qsin= aimag(c_c1_c2)   
 
724
                                   
 
725
                                   ! For COHP, insert non-trivial H matrix
 
726
                                   H_c2 = matmul( H, spinor_2 )
 
727
                                   c_c1_H_c2  = dot_product( spinor_1, H_c2 )
 
728
                                   qcos_H= real(c_c1_H_c2, dp) 
 
729
                                   qsin_H= aimag(c_c1_H_c2)
 
730
 
681
731
                                else
682
 
                                   qcos= (wf(1,io1)*wf(1,io2) + &
683
 
                                        wf(2,io1)*wf(2,io2))
684
 
                                   qsin= (wf(1,io1)*wf(2,io2) - &
685
 
                                        wf(2,io1)*wf(1,io2))
 
732
                                   ! These have explicit spin quantum numbers (is)
 
733
                                   if (gamma_wfsx) then
 
734
                                      qcos = wf(1,io1)*wf(1,io2) 
 
735
                                      qsin = 0.0_dp
 
736
                                      qcos_H = qcos * Hamilt(ind,is)
 
737
                                      qsin_H = 0.0_dp
 
738
                                   else
 
739
                                      qcos = (wf(1,io1)*wf(1,io2) + &
 
740
                                           wf(2,io1)*wf(2,io2))
 
741
                                      qcos_H = qcos * Hamilt(ind,is)
 
742
                                      qsin = (wf(1,io1)*wf(2,io2) - &
 
743
                                           wf(2,io1)*wf(1,io2))
 
744
                                      qsin_H = qsin * Hamilt(ind,is)
 
745
                                   endif
686
746
                                endif
687
747
 
688
 
                             ! k*R_12    (r_2-r_1)
689
 
                             alfa=dot_product(pk(1:3,ik),xij(1:3,ind))
 
748
                                ! k*R_12    (r_2-r_1)
 
749
                                alfa=dot_product(pk(1:3,ik),xij(1:3,ind))
690
750
 
691
 
                             ! Crb = Real(C_1*conjg(C_2)*exp(-i*alfa)) * S_12
692
 
                             !AG: This one better --  or Real(conjg(C_1)*C_2)*exp(+i*alfa)) * S_12
693
 
                             ! Common factor computed here
694
 
                             factor =  (qcos*cos(alfa)-qsin*sin(alfa)) * wk(ik)
 
751
                                ! Crb = Real(C_1*conjg(C_2)*exp(-i*alfa)) * S_12
 
752
                                ! or    Real(conjg(C_1)*C_2)*exp(+i*alfa)) * S_12
 
753
                                ! Common factor computed here
 
754
                                factor_S =  Sover(ind) * (qcos*cos(alfa)-qsin*sin(alfa)) * wk(ik)
 
755
                                factor_H =  (qcos_H*cos(alfa)-qsin_H*sin(alfa)) * wk(ik)
695
756
 
696
757
                                if (dos) then
697
758
                                   ! Note that there is no factor of 2
698
759
                                   do i = imin, imax
699
 
                                      pdos_vals(i,is,ic)=pdos_vals(i,is,ic)  +  Sover(ind)*factor*ww(i)
 
760
                                      pdos_vals(i,is,ic)=pdos_vals(i,is,ic)  +  factor_S * ww(i)
700
761
                                   enddo
701
762
                                endif
702
763
 
703
764
                                if (coop) then
704
765
                                   !! COOP is basically the off-diagonal DOS, chosen on the basis of particular bonds
705
766
                                   do i = imin, imax
706
 
                                      coop_vals(i,is,ic)=coop_vals(i,is,ic)  +  Sover(ind) * factor * ww(i)
707
 
                                      cohp_vals(i,is,ic)=cohp_vals(i,is,ic)  +  hamilt(ind,is) * factor * ww(i)
 
767
                                      coop_vals(i,is,ic)=coop_vals(i,is,ic)  +  factor_S * ww(i)
 
768
                                      cohp_vals(i,is,ic)=cohp_vals(i,is,ic)  +  factor_H * ww(i)
708
769
                                   enddo
709
770
                                endif  ! coop
710
771
 
720
781
           deallocate (list_io2)
721
782
           deallocate (list_ind)
722
783
 
723
 
     !      PDOS output
 
784
     !     Output curves
724
785
     !
725
786
           if (dos) then
726
 
              open(tab_u,file=trim(mflnm)// "." // trim(tit(ic)) // '.pdos')
727
 
              write(tab_u,"(a1,14x,'ENERGY',4(16x,a1,i1))") '#', ("s",m, m=1,nsp)
728
 
              do i=1, npts_energy
729
 
                 energy = low_e + e_step*(i-1)
730
 
                 write(tab_u,"(f20.8,10(2f13.8,5x))")  &           ! Spin in loop
731
 
                      energy, (pdos_vals(i,is,ic),is=1,nspin)
732
 
              enddo
733
 
              close(tab_u)
 
787
              call write_curve_to_file(trim(mflnm)// "." // trim(tit(ic)) // '.pdos', &
 
788
                                       pdos_vals(:,:,ic))
734
789
           endif
735
 
 
736
790
           if (coop) then
737
 
     !
738
 
     !      COOP output
739
 
     !
740
 
              open(tab_u,file=trim(mflnm)// "." // trim(tit(ic)) // '.coop')
741
 
              write(tab_u,"(a1,14x,'ENERGY',4(16x,a1,i1))") '#', ("s",m, m=1,nsp)
742
 
              do i=1, npts_energy
743
 
                 energy = low_e + e_step*(i-1)
744
 
                 write(tab_u,"(f20.8,10(2f13.8,5x))")  &           ! Spin in loop
745
 
                      energy, (coop_vals(i,is,ic),is=1,nspin)
746
 
              enddo
747
 
              close(tab_u)
748
 
              !
749
 
              !      COHP output
750
 
              !
751
 
              open(tab_u,file=trim(mflnm)// "." // trim(tit(ic)) // '.cohp')
752
 
              write(tab_u,"(a1,14x,'ENERGY',4(16x,a1,i1))") '#', ("s",m, m=1,nsp)
753
 
              do i=1, npts_energy
754
 
                 energy = low_e + e_step*(i-1)
755
 
                 write(tab_u,"(f20.8,10(2f13.8,5x))")  &           ! Spin in loop
756
 
                      energy, (cohp_vals(i,is,ic),is=1,nspin)
757
 
              enddo
758
 
              close(tab_u)
759
 
 
760
 
           endif  ! coop
761
 
 
 
791
              call write_curve_to_file(trim(mflnm)// "." // trim(tit(ic)) // '.coop', &
 
792
                                       coop_vals(:,:,ic))
 
793
              call write_curve_to_file(trim(mflnm)// "." // trim(tit(ic)) // '.cohp', &
 
794
                                       cohp_vals(:,:,ic))
 
795
           endif
 
796
           
762
797
        enddo    ! ic
763
798
 
764
799
!--------------------------------------------------------
766
801
     !
767
802
     !      Simple DOS output
768
803
     !
769
 
     call io_assign(idos)
770
 
     open(idos,file=trim(sflnm)//".ados",form="formatted", &
771
 
          status="unknown",action="write",position="rewind")
772
 
     write(idos,*) "#  Energy   SIMPLE DOS"
773
 
     e_step = (high_e-low_e)/(npts_energy-1)
774
 
     do i = 1, npts_energy
775
 
        energy = low_e + e_step*(i-1)
776
 
        ! Note division by the number of curves
777
 
        write(idos,*) energy, (ados(i,is)/ncb, is=1,nsp)
778
 
     enddo
779
 
     call io_close(idos)
780
 
 
 
804
     ! Divide by the number of curves        
 
805
     call write_curve_to_file(trim(sflnm)//".ados",ados/ncb)
781
806
 
782
807
 
783
808
CONTAINS
784
809
 
 
810
  subroutine write_curve_to_file(filename,array)
 
811
    character(len=*), intent(in) :: filename
 
812
    real(dp), intent(in) :: array(:,:)
 
813
 
 
814
    integer  :: i
 
815
    real(dp) :: energy
 
816
    
 
817
    ! tab_u, low_e, e_step, npts_energy by host association
 
818
    
 
819
    open(tab_u,file=trim(filename))
 
820
    !
 
821
    ! Header
 
822
    !
 
823
    select case ( wfs_spin_flag)
 
824
    case ( 1 )
 
825
       ! Two identical 'spin' columns, and the sum (complete value) in the third column
 
826
       write(tab_u,"(a1,12x,'ENERGY',3(14x,a))") '#', 's1', 's2=s1', 'total'
 
827
    case ( 2 )
 
828
       ! Two 'spin' columns, and the sum (complete value) in the third column
 
829
       write(tab_u,"(a1,12x,'ENERGY',3(14x,a))") '#', 's1', 's2', 'total'
 
830
    case ( 4 )
 
831
       ! A single column with  the complete value
 
832
       write(tab_u,"(a1,12x,'ENERGY',14x,a)") '#', 'total'
 
833
    end select
 
834
    !
 
835
    ! Values   
 
836
    !
 
837
    do i=1, npts_energy
 
838
       energy = low_e + e_step*(i-1)
 
839
       select case ( wfs_spin_flag)
 
840
       case ( 1 )
 
841
          ! Two identical 'spin' columns, and the sum (complete value) in the third column
 
842
          write(tab_u,"(f20.8,3(5x,f13.8))")  &   
 
843
               energy, array(i,1), array(i,1), 2*array(i,1)
 
844
       case ( 2 )
 
845
          ! Two 'spin' columns, and the sum (complete value) in the third column
 
846
          write(tab_u,"(f20.8,3(5x,f13.8))")  &      
 
847
               energy, (array(i,is),is=1,nspin_blocks), sum(array(i,:))
 
848
       case ( 4 )
 
849
          ! A single column with  the complete value
 
850
          write(tab_u,"(f20.8,5x,f13.8)") energy, array(i,1)
 
851
       end select
 
852
    enddo
 
853
    close(tab_u)
 
854
 
 
855
  end subroutine write_curve_to_file
 
856
  
785
857
  function delta(x) result(res)
786
858
    real(dp), intent(in) :: x
787
859
    real(dp)             :: res