Line data Source code
1 : !!****m* ABINIT/m_rec_tools
2 : !! NAME
3 : !! m_rec_tools
4 : !!
5 : !! FUNCTION
6 : !! This module provides some functions more or less generic used
7 : !! in the Recursion Mathod
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2002-2026 ABINIT group (MMancini)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
15 : !!
16 : !! NOTES
17 : !!
18 : !! SOURCE
19 :
20 : #if defined HAVE_CONFIG_H
21 : #include "config.h"
22 : #endif
23 :
24 : #include "abi_common.h"
25 :
26 : module m_rec_tools
27 :
28 : use defs_basis
29 : use m_abicore
30 :
31 : use defs_rectypes, only : recparall_type
32 :
33 : private
34 :
35 : public :: &
36 : get_pt0_pt1, & !--To get pt0 pt1 from inf,sup
37 : reshape_pot, & !--To rescale the potential between 2 grids
38 : trottersum !--To calculate the trotter sum
39 :
40 : CONTAINS !===========================================================
41 : !!***
42 :
43 :
44 : !!****f* m_rec_tools/get_pt0_pt1
45 : !! NAME
46 : !! get_pt0_pt1
47 : !!
48 : !! FUNCTION
49 : !! utility function to get pt0 and pt1 for given inf and sup
50 : !!
51 : !! INPUTS
52 : !! ngfft(3) = fine grid (corresponds to dtset%ngfft(1:3))
53 : !! inf = inferior point in-line coordinate on the coarse grid
54 : !! sup = superior point in-line coordinate on the coarse grid
55 : !!
56 : !! OUTPUT
57 : !! recpar%pt0<type(vec_int)>=Intial point for this proc in x,y,z
58 : !! recpar%pt1<type(vec_int)>=Final point for this proc in x,y,z
59 : !! recpar%min_pt=inferior point in-line coordinate on the fine grid
60 : !! recpar%max_pt=superior point in-line coordinate on the dine grid
61 : !!
62 : !! SIDE EFFECTS
63 : !!
64 : !! SOURCE
65 6 : subroutine get_pt0_pt1(ngfft,gratio,inf,sup,recpar)
66 :
67 : !Arguments ------------------------------------
68 : integer,intent(in) :: gratio
69 : integer,intent(in) :: inf,sup
70 : type(recparall_type),intent(inout) :: recpar
71 : integer,intent(in) :: ngfft(3)
72 : !Local ---------------------------
73 : integer :: x,y,z,count
74 : integer :: boz,boy,pt
75 : ! *********************************************************************
76 6 : count = 0
77 45 : do z = 0,ngfft(3)-1,gratio
78 45 : boz = z*ngfft(2)
79 552 : do y = 0,ngfft(2)-1,gratio
80 513 : boy = (boz+y)*ngfft(1)
81 6621 : do x = 0,ngfft(1)-1,gratio
82 6075 : pt = boy+x
83 6075 : if(count >= inf) then
84 3483 : if(count == inf) then
85 6 : recpar%pt0%x = x; recpar%pt0%y = y; recpar%pt0%z = z
86 6 : recpar%min_pt = pt
87 : end if
88 3483 : recpar%pt1%x = x; recpar%pt1%y = y; recpar%pt1%z = z
89 3483 : recpar%max_pt = pt
90 : endif
91 6075 : count = count+1
92 6582 : if(count == sup ) return
93 : end do
94 : end do
95 : end do
96 : end subroutine get_pt0_pt1
97 : !!***
98 :
99 : !!****f* m_rec_tools/reshape_pot
100 : !! NAME
101 : !! reshape_pot
102 : !!
103 : !! FUNCTION
104 : !! Reshape array on
105 : !!
106 : !! INPUTS
107 : !! nfft=size of the big grid
108 : !! nfftrec=size of the cut grid
109 : !! trasl(3) center point
110 : !! ngfft(3) dimensions of the big grid
111 : !! ngfftrec(3) dimensions of the cut grid
112 : !! pot(0:nfft-1) 3d-array on the big grid
113 : !!
114 : !! OUTPUT
115 : !! potloc is a cut of pot around trasl
116 : !!
117 : !! SOURCE
118 :
119 189 : subroutine reshape_pot(trasl,nfft,nfftrec,ngfft,ngfftrec,pot,potloc)
120 :
121 : !Arguments ------------------------------------
122 : integer, intent(in) :: nfft,nfftrec
123 : integer, intent(in) :: trasl(3)
124 : integer, intent(in) :: ngfftrec(3),ngfft(3)
125 : real(dp), intent(in) :: pot(0:nfft-1)
126 : real(dp), intent(out):: potloc(0:nfftrec-1)
127 : !Local ----------------------------------------
128 : ! scalars
129 : integer :: zz,yy,xx,parz,pary
130 : integer :: modi,modj,modk
131 : !character(len=500) :: msg
132 : ! *********************************************************************
133 945 : do zz = 0,ngfftrec(3)-1
134 756 : modk = modulo(zz-trasl(3),ngfft(3))*ngfft(2)
135 756 : parz = ngfftrec(2)*zz
136 3969 : do yy = 0,ngfftrec(2)-1
137 3024 : modj = (modulo(yy-trasl(2),ngfft(2))+modk)*ngfft(1)
138 3024 : pary = ngfftrec(1)*(yy+parz)
139 15876 : do xx = 0,ngfftrec(1)-1
140 12096 : modi = modulo(xx-trasl(1),ngfft(1))+modj
141 15120 : potloc(xx+pary) = pot(modi)
142 : end do
143 : end do
144 : end do
145 :
146 189 : end subroutine reshape_pot
147 : !!***
148 :
149 : !!****f* m_rec_tools/trottersum
150 : !! NAME
151 : !! trottersum
152 : !!
153 : !! FUNCTION
154 : !! Calculate the contribution to the partial fraction decomposition
155 : !! due to a recursion step.
156 : !!
157 : !! INPUTS
158 : !! dim_trott=dimension of the partial fraction decomposition (PFD)
159 : !! pi_on_rtrotter=parameter pi/rtrotter
160 : !! an,bn2=recursion coefficients at irec
161 : !! exp1=numerical factor (see density_rec)
162 : !! coeef_mu=numerical factor
163 : !!
164 : !! OUTPUT
165 : !!
166 : !! SIZE EFFECTS
167 : !! D,N=denominator and numerator accumalator of PFD
168 : !! Dold,Nold=denominator and numerator of PFD (old values)
169 : !! facrec0=used to select irec=0
170 : !! error=estimated error of recursion at this step
171 : !! prod_b2=numerical factor
172 : !!
173 : !! SOURCE
174 :
175 1423443 : subroutine trottersum(dim_trott,error,&
176 : & prod_b2,pi_on_rtrotter,&
177 : & facrec0,coeef_mu,exp1,&
178 : & an,bn2,&
179 1423443 : & N,D,Nold,Dold)
180 :
181 : !Arguments ------------------------------------
182 : !scalars
183 : integer, intent(in) :: dim_trott
184 : real(dp), intent(in) :: an,bn2,exp1,pi_on_rtrotter
185 : real(dp), intent(inout) :: error,prod_b2
186 : complex(dp), intent(in) :: coeef_mu
187 : complex(dp), intent(inout) :: facrec0
188 : !arrays
189 : complex(dp),intent(inout) :: D(0:dim_trott),Dold(0:dim_trott)
190 : complex(dp),intent(inout) :: N(0:dim_trott),Nold(0:dim_trott)
191 : !Local ----------------------------------------
192 : ! scalars
193 : integer :: itrot
194 : real(dp) :: arg
195 : complex(dp) :: Dnew,Nnew,zj
196 : !character(len=500) :: msg
197 : ! *********************************************************************
198 :
199 1423443 : error = zero
200 1423443 : prod_b2 = prod_b2 * exp1 * bn2
201 :
202 29892303 : do itrot=0,dim_trott
203 28468860 : arg = pi_on_rtrotter*(real( itrot,dp) + half )
204 28468860 : zj = cmplx(cos(arg),sin(arg),dp)*coeef_mu
205 :
206 : Nnew = zj*facrec0 + &
207 : & (zj - cmplx(an ,zero,dp))*N(itrot) - &
208 28468860 : & cmplx(bn2 ,zero,dp)*Nold(itrot)
209 : Dnew = (zj - cmplx(an ,zero,dp))*D(itrot) - &
210 28468860 : & cmplx(bn2 ,zero,dp)*Dold(itrot)
211 28468860 : Nold(itrot) = N(itrot)
212 28468860 : Dold(itrot) = D(itrot)
213 28468860 : N(itrot) = Nnew
214 28468860 : D(itrot) = Dnew
215 :
216 : !--Error estimator
217 29892303 : error = error + abs(prod_b2/(D(itrot)*Dold(itrot)))
218 : end do
219 1423443 : facrec0 = czero
220 :
221 1423443 : end subroutine trottersum
222 : !!***
223 :
224 : end module m_rec_tools
225 : !!***
|
