Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2023 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_dfpt_hsvac
7 : USE m_juDFT
8 : CONTAINS
9 : !-----------------------------------------------------------------------------
10 : ! Calculate the vacuum contribution to the Hamiltonian perturbation
11 : !-----------------------------------------------------------------------------
12 0 : SUBROUTINE dfpt_hsvac(vacuum, stars, fmpi, jsp, input, v, v1, evac, cell, &
13 0 : & lapwq, lapw, noco, nococonv, hmat)
14 :
15 : USE m_vacfun
16 : USE m_types
17 :
18 : IMPLICIT NONE
19 :
20 : TYPE(t_input),INTENT(IN) :: input
21 : TYPE(t_vacuum),INTENT(IN) :: vacuum
22 : TYPE(t_noco),INTENT(IN) :: noco
23 : TYPE(t_nococonv),INTENT(IN) :: nococonv
24 : TYPE(t_stars),INTENT(IN) :: stars
25 : TYPE(t_cell),INTENT(IN) :: cell
26 : TYPE(t_lapw),INTENT(IN) :: lapw, lapwq
27 : TYPE(t_mpi),INTENT(IN) :: fmpi
28 : TYPE(t_potden),INTENT(IN) :: v, v1
29 : CLASS(t_mat),INTENT(INOUT) :: hmat(:,:)
30 :
31 : INTEGER, INTENT (IN) :: jsp
32 : REAL, INTENT (IN) :: evac(2,input%jspins)
33 :
34 : COMPLEX :: hij,c_1
35 : REAL :: d2,gz,sign,th,wronk
36 : INTEGER :: ikG,ikG2,ikGPr,ikGPr2,jspin,ikG0
37 : INTEGER :: ivac,igSpin,igSpinPr
38 : INTEGER :: iSpin,iSpinPr
39 :
40 0 : INTEGER :: nv2(input%jspins)
41 0 : INTEGER :: kvac(2,lapw%dim_nv2d(),input%jspins)
42 0 : INTEGER :: map2(lapw%dim_nvd(),input%jspins)
43 0 : COMPLEX :: tddv(lapw%dim_nv2d(),lapw%dim_nv2d()),tduv(lapw%dim_nv2d(),lapw%dim_nv2d())
44 0 : COMPLEX :: tudv(lapw%dim_nv2d(),lapw%dim_nv2d()),tuuv(lapw%dim_nv2d(),lapw%dim_nv2d())
45 0 : COMPLEX :: a(lapw%dim_nvd(),input%jspins),b(lapw%dim_nvd(),input%jspins)
46 0 : REAL :: ddnv(lapw%dim_nv2d(),input%jspins),dudz(lapw%dim_nv2d(),input%jspins)
47 0 : REAL :: duz(lapw%dim_nv2d(),input%jspins), udz(lapw%dim_nv2d(),input%jspins)
48 0 : REAL :: uz(lapw%dim_nv2d(),input%jspins), dudzq(lapw%dim_nv2d(),input%jspins)
49 0 : REAL :: duzq(lapw%dim_nv2d(),input%jspins), udzq(lapw%dim_nv2d(),input%jspins)
50 0 : REAL :: uzq(lapw%dim_nv2d(),input%jspins)
51 0 : COMPLEX :: aq(lapw%dim_nvd(),input%jspins),bq(lapw%dim_nvd(),input%jspins)
52 0 : INTEGER :: map2q(lapw%dim_nvd(),input%jspins)
53 0 : INTEGER :: nv2q(input%jspins)
54 0 : INTEGER :: kvacq(2,lapw%dim_nv2d(),input%jspins)
55 :
56 0 : d2 = SQRT(cell%omtil/cell%area)
57 :
58 : !---> set up mapping function from 3d-->2d lapws
59 0 : DO jspin = 1,input%jspins
60 0 : nv2(jspin) = 0
61 0 : k_loop:DO ikG = 1, lapw%nv(jspin)
62 0 : DO ikG2 = 1, nv2(jspin)
63 0 : IF (all(lapw%gvec(1:2,ikG,jspin)==kvac(1:2,ikG2,jspin))) THEN
64 0 : map2(ikG,jspin) = ikG2
65 0 : CYCLE k_loop
66 : END IF
67 : END DO
68 0 : nv2(jspin) = nv2(jspin) + 1
69 0 : IF (nv2(jspin)>lapw%dim_nv2d()) CALL juDFT_error("hsvac:lapw%dim_nv2d()",calledby ="hsvac")
70 0 : kvac(1:2,nv2(jspin),jspin) = lapw%gvec(1:2,ikG,jspin)
71 0 : map2(ikG,jspin) = nv2(jspin)
72 : END DO k_loop
73 : END DO
74 :
75 0 : DO jspin = 1,input%jspins
76 0 : nv2q(jspin) = 0
77 0 : k_loopq:DO ikG = 1, lapwq%nv(jspin)
78 0 : DO ikG2 = 1, nv2q(jspin)
79 0 : IF (all(lapwq%gvec(1:2,ikG,jspin)==kvacq(1:2,ikG2,jspin))) THEN
80 0 : map2q(ikG,jspin) = ikG2
81 0 : CYCLE k_loopq
82 : END IF
83 : END DO
84 0 : nv2q(jspin) = nv2q(jspin) + 1
85 0 : IF (nv2q(jspin)>lapw%dim_nv2d()) CALL juDFT_error("hsvac:lapw%dim_nv2d()",calledby ="hsvac")
86 0 : kvacq(1:2,nv2q(jspin),jspin) = lapwq%gvec(1:2,ikG,jspin)
87 0 : map2q(ikG,jspin) = nv2q(jspin)
88 : END DO k_loopq
89 : END DO
90 :
91 : !---> loop over the two vacuua (1: upper; 2: lower)
92 0 : DO ivac = 1,2
93 0 : sign = 3. - 2.*ivac !+/- 1
94 0 : DO iSpin=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin
95 0 : igSpin=MIN(SIZE(hmat,1),iSpin) !in colinear case igSpin=1
96 0 : DO iSpinPr=MERGE(1,jsp,noco%l_noco),MERGE(2,jsp,noco%l_noco) !loop over global spin
97 0 : igSpinPr=MIN(SIZE(hmat,1),iSpinPr) !in colinear case igSpinPr=1
98 : !---> get the wavefunctions and set up the tuuv, etc matrices
99 0 : CALL timestart("vacfun")
100 : CALL vacfun(fmpi, vacuum, stars, input, nococonv, iSpin, iSpinPr, &
101 : & cell, ivac, evac, lapw%bkpt, v%vac(:vacuum%nmzxyd,2:,:,:), v%vac(:,1,:,:), kvac, nv2, &
102 : & tuuv, tddv, tudv, tduv, uz, duz, udz, dudz, ddnv, wronk,&
103 0 : & lapwq%bkpt+lapwq%qphon, v1%vac, kvacq, nv2q, uzq, duzq, udzq, dudzq)
104 0 : CALL timestop("vacfun")
105 :
106 : !---> generate a and b coeffficients
107 0 : DO jspin = MIN(iSpin,iSpinPr),MAX(iSpin,iSpinPr)
108 0 : DO ikG = 1,lapw%nv(jspin)
109 0 : gz = sign*cell%bmat(3,3)*lapw%k3(ikG,jspin)
110 0 : ikG2 = map2(ikG,jspin)
111 0 : th = gz*cell%z1
112 0 : c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
113 0 : a(ikG,jspin) = - c_1 * CMPLX(dudz(ikG2,jspin), gz*udz(ikG2,jspin) )
114 0 : b(ikG,jspin) = c_1 * CMPLX(duz(ikG2,jspin), gz* uz(ikG2,jspin) )
115 : END DO
116 :
117 0 : DO ikG = 1,lapwq%nv(jspin)
118 0 : gz = sign*cell%bmat(3,3)*lapwq%k3(ikG,jspin)
119 0 : ikG2 = map2q(ikG,jspin)
120 0 : th = gz*cell%z1
121 0 : c_1 = CMPLX( COS(th), SIN(th) )/ (d2*wronk)
122 0 : aq(ikG,jspin) = - c_1 * CMPLX(dudzq(ikG2,jspin), gz*udzq(ikG2,jspin) )
123 0 : bq(ikG,jspin) = c_1 * CMPLX(duzq(ikG2,jspin), gz* uzq(ikG2,jspin) )
124 : END DO
125 : END DO
126 :
127 : !---> hamiltonian update
128 0 : DO ikG = fmpi%n_rank+1,lapw%nv(iSpin),fmpi%n_size
129 0 : ikG0 = (ikG-1)/fmpi%n_size + 1 !local column index
130 0 : ikG2 = map2(ikG,iSpin)
131 0 : DO ikGPr = 1, lapwq%nv(iSpinPr)
132 0 : ikGPr2 = map2q(ikGPr, iSpinPr)
133 : hij = CONJG(aq(ikGPr, iSpinPr)) * tuuv(ikGPr2, ikG2) * a(ikG,iSpin) &
134 : + CONJG(bq(ikGPr, iSpinPr)) * tddv(ikGPr2, ikG2) * b(ikG,iSpin) &
135 : + CONJG(aq(ikGPr, iSpinPr)) * tudv(ikGPr2, ikG2) * b(ikG,iSpin) &
136 0 : + CONJG(bq(ikGPr, iSpinPr)) * tduv(ikGPr2, ikG2) * a(ikG,iSpin)
137 0 : hmat(igSpinPr,igSpin)%data_c(ikGPr,ikG0) = hmat(igSpinPr,igSpin)%data_c(ikGPr,ikG0) + hij
138 : END DO
139 : END DO
140 : !---> end of loop over different parts of the potential matrix
141 : END DO
142 : !---> end of loop over vacua
143 : END DO
144 : END DO
145 0 : END SUBROUTINE dfpt_hsvac
146 : END MODULE m_dfpt_hsvac
|
