Branch data 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 : : !!***
|
