TD-RAS-SCF/module_index_radial.f90 at master · danielreich/TD-RAS-SCF · GitHub

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

!
! the spare operator are using in the calcuations
!
 module not_full_operator_index
  use global
  use fedvr3d_basis_set
  implicit none
  integer :: n_total_kinetic  ! the number of kinetic element with non-zero
  integer,allocatable,save :: index_kinetic_basis(:,:)
  integer,allocatable,save :: index_ham_act(:,:)
  integer,allocatable,save :: index_kinetic_operator(:,:)
  contains


  subroutine index_of_kinetic
   implicit none
   integer :: idicp,isum,ibasis,jbasis,e,f,itemp,jtemp
   integer :: i,j
!---------------------------------------------------------------
!                radial part
! figure out how many element in kinetic operator are non-zero
!---------------------------------------------------------------
   isum = 0
   isum = (fedvr3d%fedvr_nb(1) -1)**2

   do idicp =2, fedvr3d%number_of_element -1
     isum = isum + (fedvr3d%fedvr_nb(idicp))**2
   enddo

    isum = isum + (fedvr3d%fedvr_nb(fedvr3d%number_of_element) -1)**2

    n_total_kinetic = isum - (fedvr3d%number_of_element - 1)

    n_total_kinetic = (n_total_kinetic - fedvr3d%nb_r)/2 + fedvr3d%nb_r
!
! there are total n_total_kinetic non-zero elements in radial part kinetic operator matrix element
!
   allocate( index_kinetic_basis(n_total_kinetic,2)  )
   allocate( index_ham_act(fedvr3d%nb_r,2))
   allocate(index_kinetic_operator(fedvr3d%nb_r,fedvr3d%nb_r))

    itemp = 0
    do ibasis = 1, fedvr3d%nb_r

     ! jtemp = 0
      do jbasis =1,ibasis

       e = which_element(ibasis) ; i = which_basis(ibasis)
       f = which_element(jbasis) ; j = which_basis(jbasis)


      if(e==f .or. ((e==(f-1)).and. (i==fedvr3d%fedvr_nb(e))) .or.((e==(f+1)).and.(j==fedvr3d%fedvr_nb(f)))) then
        itemp = itemp +1
        index_kinetic_basis(itemp,1) = ibasis
        index_kinetic_basis(itemp,2) = jbasis

        index_kinetic_operator(ibasis,jbasis) = itemp


      !  jtemp = jtemp +1
      !  if(jtemp ==1) then
      !      index_ham_act(ibasis,1) = jbasis
      !  endif
      !  index_ham_act(ibasis,2) = jbasis
       endif
      enddo
   enddo

!
!=================================================================================
!
    do ibasis = 1, fedvr3d%nb_r

      jtemp = 0
      do jbasis =1,fedvr3d%nb_r

       e = which_element(ibasis) ; i = which_basis(ibasis)
       f = which_element(jbasis) ; j = which_basis(jbasis)
      if(e==f .or. ((e==(f-1)).and. (i==fedvr3d%fedvr_nb(e))) .or.((e==(f+1)).and.(j==fedvr3d%fedvr_nb(f)))) then
        jtemp = jtemp +1
if(jtemp ==1) then
  index_ham_act(ibasis,1) = jbasis
endif
          index_ham_act(ibasis,2) = jbasis
      endif
      enddo
   enddo


    if(itemp /= n_total_kinetic) then
       write(*,*) 'error happen in set_kinetic_energy_operator'
       stop
    endif

   return
  end subroutine index_of_kinetic


 end module not_full_operator_index

!
! the full
!
subroutine drive_not_full_operator_index
 use not_full_operator_index
 implicit none

 call index_of_kinetic

return
end subroutine drive_not_full_operator_index