← Back to branch summary

~albertog/siesta/efield-2.5-jjunquera

« back to all changes in this revision

Viewing changes to Util/COOP/orbital_set.f90

  • Committer: Alberto Garcia
  • Date: 2007-11-22 21:31:50 UTC
  • mfrom: (unknown (missing))
  • Revision ID: Arch-1:siesta@uam.es--2006%siesta-devel--reference--2.5--patch-1
Merged COOP/COHP functionality
Merged the COOP/COHP/PDOS code from the 2.1 branch complex.

��-- Functionality changes:

* The wavefunction file is now in WFSX format (well-packed, single-precision, and thus
significantly smaller). The Util/wfsx2wfs program can be used to re-generate the old format.

* The xijo array (for relative positions between interacting orbitals) is always produced.

-- New features:

Option COOP.Write triggers the creation of a wavefunction file (for the SCF-k-point set)
and an HSX file (enhanced HS format) that can later be processed by Util/COOP/mprop
for the off-line generation of COOP/COHP/(P)DOS.

(replayed patch-1 of outlier branch ref-2.4)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Util/COOP/orbital_set.f90
 
1
module orbital_set
 
2
CONTAINS
 
3
subroutine get_orbital_set(line,set_mask)
 
4
  use main_vars
 
5
  use subs, only : txt2wrd, orbital
 
6
 
 
7
  implicit none
 
8
 
 
9
  character(len=*), intent(in) :: line
 
10
  logical, intent(out)         :: set_mask(:)
 
11
 
 
12
  print *, "Size of set_mask: ", size(set_mask)
 
13
  set_mask(:) = .false.
 
14
 
 
15
  call txt2wrd (line, wrd, nwd, nlwmx)
 
16
  if (nwd.gt.nlwmx) stop "* Groups per subset limit exceeded."
 
17
 
 
18
  if (trim(wrd(1)).eq.'+') then
 
19
     do iw=2,nwd
 
20
        k=ival(wrd(iw))
 
21
        if (k.le.0.or.k.gt.no_u) then
 
22
           print *, "Wrong orbital number: ", k
 
23
           STOP
 
24
        endif
 
25
        set_mask(k) = .true.
 
26
     enddo
 
27
 
 
28
  else
 
29
 
 
30
     do iw=1,nwd
 
31
        call orbital (wrd(iw), ia, cx, n, l, k)
 
32
 
 
33
        if (ia.lt.0) then
 
34
           print *, "Error in orb spec: ", trim(wrd(iw))
 
35
           STOP
 
36
        endif
 
37
 
 
38
        if (ia.eq.0) then
 
39
           it=0
 
40
           do i=1,nspecies
 
41
              if (trim(label(i)) .eq. trim(cx)) it=i
 
42
           enddo
 
43
           if (it.eq.0) then
 
44
              print *, "Wrong species: ", trim(cx)
 
45
              STOP
 
46
           endif
 
47
        endif
 
48
        if (ia > na_u) then
 
49
           print *, "Atom index too big: ", ia
 
50
           STOP
 
51
        endif
 
52
        !!! Could go on checking whether a given
 
53
        !!! atom has the orbitals specified, instead
 
54
        !!! or giving an empty result
 
55
        !
 
56
        !             See which orbitals match
 
57
        !
 
58
        do io=1,no_u
 
59
           if ((za(io).eq.ia).or.(ia.eq.0.and.zc(io).eq.it)) then
 
60
              if ((zn(io).eq.n.or.n.eq.-1).and. &
 
61
                   (zl(io).eq.l.or.l.eq.-1).and. &
 
62
                   (zx(io).eq.k.or.k.eq.-1)) then
 
63
                 set_mask(io) = .true.
 
64
              endif
 
65
           endif
 
66
        enddo
 
67
     enddo
 
68
 
 
69
  endif
 
70
 
 
71
end subroutine get_orbital_set
 
72
 
 
73
end module orbital_set
 
74