Line data Source code
1 : !!****m* ABINIT/m_macroave
2 : !! NAME
3 : !! m_macroave
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1997-2003 J.Soler, P. Ordejon, J. Junquera
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_macroave
22 :
23 : implicit none
24 :
25 : ! private
26 :
27 : public :: POLINT
28 : !!***
29 :
30 : contains
31 : !!***
32 :
33 0 : subroutine iorho( task, fname, cell, mesh, nsm, maxp, nspin, f, found )
34 :
35 : ! *********************************************************************
36 : ! Saves/recovers the electron density at the mesh points.
37 : ! This simplified version reads only files written in serial mode.
38 :
39 : ! Writen by J.Soler July 1997.
40 : ! Copyright (C) 1997-2003 J.Soler, P. Ordejon, J. Junquera
41 : ! This file is distributed under the terms of the
42 : ! GNU General Public License, see ~abinit/COPYING
43 : ! or http://www.gnu.org/copyleft/gpl.txt .
44 :
45 : ! *************************** INPUT **********************************
46 : ! character*(*) task : 'read'/'READ' or 'write'/'WRITE'
47 : ! character*(*) fname : File name for input or output
48 : ! integer nsm : Number of sub-mesh points per mesh point
49 : ! (not used in this version)
50 : ! integer maxp : First dimension of array rho
51 : ! integer nspin : Second dimension of array rho
52 : ! ************************** OUTPUT **********************************
53 : ! integer maxp : Required first dimension of array rho,
54 : ! equal to mesh(1)*mesh(2)*mesh(3)
55 : ! Set only when task='read' and required
56 : ! value is larger than input value
57 : ! integer nspin : Number of spin polarizations (1 or 2)
58 : ! logical found : Were data found? (only when task='read')
59 : ! ******************** INPUT or OUTPUT (depending on task) ***********
60 : ! real*8 cell(3,3) : Lattice vectors
61 : ! integer mesh(3) : Number of mesh divisions of each
62 : ! lattice vector
63 : ! real f(maxp,nspin) : Electron density
64 : ! Notice single precision in this version
65 : ! *************************** UNITS ***********************************
66 : ! Units should be consistent between task='read' and 'write'
67 : ! ******************** BEHAVIOUR **************************************
68 : ! If task='read', and the values of maxp or nspin on input are less than
69 : ! those required to copy the array f from the file, then the required
70 : ! values of maxp and nspin are returned on output, but f is not read.
71 : ! *********************************************************************
72 :
73 : ! Arguments
74 : character*(*) fname, task
75 : integer maxp, mesh(3), nspin, nsm
76 : real f(maxp,nspin)
77 : real(kind=kind(0.0d0)) cell(3,3)
78 : logical found
79 :
80 : ! Internal variables and arrays
81 : character(len=11) :: fform
82 : integer i2, i3, ind, ip, is, np, ns
83 :
84 : if(.false.)write(11,*)nsm
85 : if(.false.)write(11,*)task
86 :
87 : ! Fix whether formatted or unformatted files will be used
88 0 : fform = 'unformatted'
89 :
90 : ! Look for data file
91 0 : inquire( file=fname, exist=found )
92 0 : if (.not.found) return
93 :
94 : ! Read unit cell vectors, number of mesh points and spin components
95 0 : open( unit=1, file=fname, status='old', form=fform )
96 0 : if (fform == 'formatted') then
97 0 : read(1,*) cell
98 0 : read(1,*) mesh, ns
99 : else
100 0 : read(1,*) cell
101 0 : read(1,*) mesh, ns
102 : endif
103 :
104 : ! Read density (only if array f is large enough)
105 0 : np = mesh(1) * mesh(2) * mesh(3)
106 0 : if (ns>nspin .or. np>maxp) then
107 0 : maxp = np
108 : else
109 0 : if (fform == 'formatted') then
110 : ind = 0
111 0 : do is = 1,ns
112 0 : do i3 = 1,mesh(3)
113 0 : do i2 = 1,mesh(2)
114 0 : read(1,*) (f(ind+ip,is),ip=1,mesh(1))
115 0 : ind = ind + mesh(1)
116 : enddo
117 : enddo
118 : enddo
119 : else
120 : ind = 0
121 0 : do is = 1,ns
122 0 : do i3 = 1,mesh(3)
123 0 : do i2 = 1,mesh(2)
124 0 : read(1) (f(ind+ip,is),ip=1,mesh(1))
125 0 : ind = ind + mesh(1)
126 : enddo
127 : enddo
128 : enddo
129 : endif
130 : endif
131 0 : close(1)
132 0 : nspin = ns
133 : end SUBROUTINE iorho
134 : !!***
135 :
136 4 : SUBROUTINE FOUR1(DATA,NN,ISIGN)
137 : !**********************************************************************
138 : ! Discrete Fourier transform.
139 : !**********************************************************************
140 : ! Input:
141 : ! real*8 DATA(2*NN) : Function to be Fourier transformed
142 : ! integer NN : Number of points. Must be a power of 2
143 : ! integer ISIGN : ISIG=+1/-1 => Direct/inverse transform
144 : ! Output:
145 : ! real*8 DATA(2*NN) : Fourier transformed function
146 : !**********************************************************************
147 : INTEGER :: NN, ISIGN
148 : real(kind=kind(0.0d0)) :: DATA(2*NN)
149 :
150 : INTEGER :: I, ISTEP, J, M, MMAX, N
151 : real(kind=kind(0.0d0)) :: TEMPI, TEMPR, THETA, WI, WPI, WPR, WR, WTEMP
152 : DOUBLE PRECISION, PARAMETER :: TWOPI=6.28318530717959D0,&
153 : & HALF=0.5D0, ONE=1.D0, TWO=2.D0, ZERO=0.D0
154 :
155 4 : N=2*NN
156 4 : J=1
157 :
158 4 : DO I=1,N,2
159 16384 : IF(J>I)THEN
160 8064 : TEMPR=DATA(J)
161 8064 : TEMPI=DATA(J+1)
162 8064 : DATA(J)=DATA(I)
163 8064 : DATA(J+1)=DATA(I+1)
164 8064 : DATA(I)=TEMPR
165 8064 : DATA(I+1)=TEMPI
166 : ENDIF
167 16384 : M=N/2
168 16380 : DO ! until following condition is met
169 32764 : IF ((M<2).OR.(J<=M)) EXIT
170 16380 : J=J-M
171 16380 : M=M/2
172 : END DO
173 16384 : J=J+M
174 : END DO ! I
175 : MMAX=2
176 : DO ! until following condition is met
177 52 : IF (N<=MMAX) EXIT
178 48 : ISTEP=2*MMAX
179 48 : THETA=TWOPI/(ISIGN*MMAX)
180 48 : WPR=(-TWO)*SIN(HALF*THETA)**2
181 48 : WPI=SIN(THETA)
182 48 : WR=ONE
183 48 : WI=ZERO
184 48 : DO M=1,MMAX,2
185 16380 : DO I=M,N,ISTEP
186 98304 : J=I+MMAX
187 98304 : TEMPR=WR*DATA(J)-WI*DATA(J+1)
188 98304 : TEMPI=WR*DATA(J+1)+WI*DATA(J)
189 98304 : DATA(J)=DATA(I)-TEMPR
190 98304 : DATA(J+1)=DATA(I+1)-TEMPI
191 98304 : DATA(I)=DATA(I)+TEMPR
192 98304 : DATA(I+1)=DATA(I+1)+TEMPI
193 : END DO ! I
194 16380 : WTEMP=WR
195 16380 : WR=WR*WPR-WI*WPI+WR
196 16380 : WI=WI*WPR+WTEMP*WPI+WI
197 : END DO ! M
198 4 : MMAX=ISTEP
199 : END DO ! until (N<=MMAX)
200 :
201 4 : END SUBROUTINE FOUR1
202 : !!***
203 :
204 0 : SUBROUTINE POLINT(XA,YA,N,X,Y,DY)
205 : !*****************************************************************
206 : ! Polinomic interpolation.
207 : ! D. Sanchez-Portal, Oct. 1996
208 : !*****************************************************************
209 : ! Input:
210 : ! real*8 XA(N) : x values of the function y(x) to interpolate
211 : ! real*8 YA(N) : y values of the function y(x) to interpolate
212 : ! integer N : Number of data points
213 : ! real*8 X : x value at which the interpolation is desired
214 : ! Output:
215 : ! real*8 Y : interpolated value of y(x) at X
216 : ! real*8 DY : accuracy estimate
217 : !*****************************************************************
218 : INTEGER :: N
219 : real(kind=kind(0.0d0)) :: XA(N),YA(N), X, Y, DY
220 :
221 : INTEGER :: I, M, NS
222 0 : real(kind=kind(0.0d0)) :: C(N), D(N), DEN, DIF, DIFT, HO, HP, W
223 : DOUBLE PRECISION, PARAMETER :: ZERO=0.D0
224 :
225 0 : NS=1
226 0 : DIF=ABS(X-XA(1))
227 0 : DO I=1,N
228 0 : DIFT=ABS(X-XA(I))
229 0 : IF (DIFT<DIF) THEN
230 0 : NS=I
231 0 : DIF=DIFT
232 : ENDIF
233 0 : C(I)=YA(I)
234 0 : D(I)=YA(I)
235 : END DO ! I
236 0 : Y=YA(NS)
237 0 : NS=NS-1
238 0 : DO M=1,N-1
239 0 : DO I=1,N-M
240 0 : HO=XA(I)-X
241 0 : HP=XA(I+M)-X
242 0 : W=C(I+1)-D(I)
243 0 : DEN=HO-HP
244 0 : IF (DEN==ZERO) STOP 'polint: ERROR. Two XAs are equal'
245 0 : DEN=W/DEN
246 0 : D(I)=HP*DEN
247 0 : C(I)=HO*DEN
248 : END DO ! I
249 0 : IF (2*NS<N-M) THEN
250 0 : DY=C(NS+1)
251 : ELSE
252 0 : DY=D(NS)
253 0 : NS=NS-1
254 : ENDIF
255 0 : Y=Y+DY
256 : END DO ! M
257 :
258 0 : END SUBROUTINE POLINT
259 :
260 3 : SUBROUTINE MACROAV_SPLINE(DX,Y,N,YP1,YPN,Y2)
261 : !***********************************************************
262 : ! Cubic Spline Interpolation.
263 : ! D. Sanchez-Portal, Oct. 1996.
264 : ! Input:
265 : ! real*8 DX : x interval between data points
266 : ! real*8 Y(N) : value of y(x) at data points
267 : ! integer N : number of data points
268 : ! real*8 YP1 : value of dy/dx at X1 (first point)
269 : ! real*8 YPN : value of dy/dx at XN (last point)
270 : ! Output:
271 : ! real*8 Y2(N): array to be used by routine MACROAV_SPLINT
272 : ! Behavior:
273 : ! - If YP1 or YPN are larger than 1E30, the natural spline
274 : ! condition (d2y/dx2=0) at the corresponding edge point.
275 : !************************************************************
276 : INTEGER :: N
277 : real(kind=kind(0.0d0)) :: DX, Y(N), YP1, YPN, Y2(N)
278 :
279 : INTEGER :: I, K
280 6 : real(kind=kind(0.0d0)) :: QN, P, SIG, U(N), UN
281 : DOUBLE PRECISION, PARAMETER :: YPMAX=0.99D30, &
282 : & HALF=0.5D0, ONE=1.D0, THREE=3.D0, TWO=2.D0, ZERO=0.D0
283 :
284 3 : IF (YP1>YPMAX) THEN
285 0 : Y2(1)=ZERO
286 0 : U(1)=ZERO
287 : ELSE
288 3 : Y2(1)=-HALF
289 3 : U(1)=(THREE/DX)*((Y(2)-Y(1))/DX-YP1)
290 : ENDIF
291 160 : DO I=2,N-1
292 157 : SIG=HALF
293 157 : P=SIG*Y2(I-1)+TWO
294 157 : Y2(I)=(SIG-ONE)/P
295 : U(I)=(THREE*( Y(I+1)+Y(I-1)-TWO*Y(I) )/(DX*DX)&
296 160 : & -SIG*U(I-1))/P
297 : END DO ! I
298 3 : IF (YPN>YPMAX) THEN
299 : QN=ZERO
300 : UN=ZERO
301 : ELSE
302 3 : QN=HALF
303 3 : UN=(THREE/DX)*(YPN-(Y(N)-Y(N-1))/DX)
304 : ENDIF
305 3 : Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+ONE)
306 163 : DO K=N-1,1,-1
307 163 : Y2(K)=Y2(K)*Y2(K+1)+U(K)
308 : END DO ! K
309 :
310 3 : END SUBROUTINE MACROAV_SPLINE
311 :
312 :
313 12288 : SUBROUTINE MACROAV_SPLINT(DX,YA,Y2A,N,X,Y,DYDX)
314 : !***************************************************************
315 : ! Cubic Spline Interpolation.
316 : ! D. Sanchez-Portal, Oct. 1996.
317 : ! Input:
318 : ! real*8 DX : x interval between data points
319 : ! real*8 YA(N) : value of y(x) at data points
320 : ! real*8 Y2A(N): array returned by routine MACROAV_SPLINE
321 : ! integer N : number of data points
322 : ! real*8 X : point at which interpolation is desired
323 : ! real*8 Y : interpolated value of y(x) at point X
324 : ! real*8 DYDX : interpolated value of dy/dx at point X
325 : !***************************************************************
326 : INTEGER :: N
327 : real(kind=kind(0.0d0)) :: DX, YA(N), Y2A(N), X, Y, DYDX
328 :
329 : INTEGER :: NHI, NLO
330 : real(kind=kind(0.0d0)) :: A, B
331 : DOUBLE PRECISION, PARAMETER ::&
332 : & ONE=1.D0, THREE=3.D0, SIX=6.D0, ZERO=0.D0
333 :
334 12288 : IF (DX==ZERO) STOP 'splint: ERROR: DX=0'
335 12288 : NLO=INT(X/DX)+1
336 12288 : NHI=NLO+1
337 12288 : A=NHI-X/DX-1
338 12288 : B=ONE-A
339 : Y=A*YA(NLO)+B*YA(NHI)+&
340 12288 : & ((A**3-A)*Y2A(NLO)+(B**3-B)*Y2A(NHI))*(DX**2)/SIX
341 : DYDX=(YA(NHI)-YA(NLO))/DX+&
342 : & (-((THREE*(A**2)-ONE)*Y2A(NLO))+&
343 12288 : & (THREE*(B**2)-ONE)*Y2A(NHI))*DX/SIX
344 :
345 12288 : END SUBROUTINE MACROAV_SPLINT
346 :
347 : ! Copyright (C) 1999-2003 (P. Ordejon, J. Junquera)
348 : ! This file is distributed under the terms of the
349 : ! GNU General Public License, see ~abinit/COPYING
350 : ! or http://www.gnu.org/copyleft/gpl.txt .
351 :
352 4 : real(kind=kind(0.0d0)) FUNCTION SURPLA( C )
353 :
354 : ! CALCULATES THE SRFACE OF THE UNIT CELL NORMAL TO THE INTERFACE
355 :
356 : real(kind=kind(0.0d0)) C(3,3)
357 : SURPLA = ( C(2,1)*C(3,2) - C(3,1)*C(2,2) ) **2 +&
358 : & ( C(3,1)*C(1,2) - C(1,1)*C(3,2) ) **2 +&
359 4 : & ( C(1,1)*C(2,2) - C(2,1)*C(1,2) ) **2
360 4 : SURPLA = SQRT( ABS( SURPLA ) )
361 4 : END FUNCTION SURPLA
362 :
363 5 : subroutine thetaft(n,L,lav,ft)
364 :
365 : ! Copyright (C) 1999-2003 (P. Ordejon, J. Junquera)
366 : ! This file is distributed under the terms of the
367 : ! GNU General Public License, see ~abinit/COPYING
368 : ! or http://www.gnu.org/copyleft/gpl.txt .
369 :
370 : !
371 : ! Calculates the Fourier coefficients of the convoluting function
372 : ! w(x) = (1/lav) theta(lav/2 - abs(x))
373 : ! asuming it is periodic with period L:
374 : !
375 : ! | |
376 : ! 1/lav|__________ __________|
377 : ! | | | |
378 : ! |__________|_____________________|_________|___
379 : ! 0 lav/2 L-lav/2 L
380 : !
381 : integer n,j
382 : real(kind=kind(0.0d0)) :: L,lav
383 : real(kind=kind(0.0d0)) :: ft(2*n)
384 : real(kind=kind(0.0d0)) :: pi
385 :
386 5 : pi=4.0d0*datan(1.0d0)
387 :
388 5 : ft(1)=n/L
389 5 : ft(2)=0.0
390 :
391 10245 : do j=2,n/2+1
392 : ft(2*(j-1)+1) = (n/(pi*lav*(j-1)/2.))*&
393 : & sin(pi*lav*(j-1)/2./L)*&
394 : & cos(pi*(L-lav/2.)*(j-1)/L)*&
395 10240 : & cos(pi*L*(j-1)/L)
396 : ft(2*(j-1)+2) = (n/(pi*lav*(j-1)/2.))*&
397 : & sin(pi*lav*(j-1)/2./L)*&
398 : & cos(pi*(L-lav/2.)*(j-1)/L)*&
399 10245 : & sin(pi*L*(j-1)/L)
400 : enddo
401 :
402 10240 : do j=n/2+2,n
403 : ft(2*(j-1)+1) = (n/(pi*lav*(-n+j-1)/2.))*&
404 : & sin(pi*lav*(-n+j-1)/2./L)*&
405 : & cos(pi*(L-lav/2.)*(-n+j-1)/L)*&
406 10235 : & cos(pi*L*(-n+j-1)/L)
407 : ft(2*(j-1)+2) = (n/(pi*lav*(-n+j-1)/2.))*&
408 : & sin(pi*lav*(-n+j-1)/2./L)*&
409 : & cos(pi*(L-lav/2.)*(-n+j-1)/L)*&
410 10240 : & sin(pi*L*(-n+j-1)/L)
411 : enddo
412 :
413 :
414 5 : return
415 : end subroutine thetaft
416 :
417 4 : DOUBLE PRECISION FUNCTION VOLCEL( C )
418 :
419 : ! Copyright (C) 1999-2003 (P. Ordejon, J. Junquera)
420 : ! This file is distributed under the terms of the
421 : ! GNU General Public License, see ~abinit/COPYING
422 : ! or http://www.gnu.org/copyleft/gpl.txt .
423 :
424 : ! CALCULATES THE VOLUME OF THE UNIT CELL
425 : DOUBLE PRECISION C(3,3)
426 : VOLCEL = ( C(2,1)*C(3,2) - C(3,1)*C(2,2) ) * C(1,3) +&
427 : & ( C(3,1)*C(1,2) - C(1,1)*C(3,2) ) * C(2,3) +&
428 4 : & ( C(1,1)*C(2,2) - C(2,1)*C(1,2) ) * C(3,3)
429 4 : VOLCEL = ABS( VOLCEL )
430 4 : END FUNCTION VOLCEL
431 :
432 : end module m_macroave
433 : !!***
|
