Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 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_grdrsis
7 : USE m_juDFT
8 : PRIVATE
9 : INTEGER,PARAMETER :: ndvgrd=6 ! this should be consistent across GGA derivative routines
10 : PUBLIC grdrsis
11 : !.....-----------------------------------------------------------------
12 : ! evaluates gradient of an interstitial function in real space
13 : !
14 : ! based on 'rhzgrd' coded by t.asada. june,1995.
15 : !.....-----------------------------------------------------------------
16 : CONTAINS
17 496 : SUBROUTINE grdrsis(ro,cell,xmax1,xmax2,xmax3, dro)
18 : !.....-----------------------------------------------------------------
19 : ! input:
20 : ! ro(0:xmax1*xmax2*xmax3-1)
21 : ! any quantity stored in usual interst. box (xmax1 x xmax2 x xmax3)
22 : ! bmat
23 : ! bravais matrix of reciprocal space
24 : ! ndvgrd
25 : ! number of ponts used when calculating derivative (3 <= ndvgrd <= 6)
26 : !
27 : ! output:
28 : ! dro(0:xmax1*xmax2*xmax3-1,3)
29 : ! gradient of ro in non-internal coordinates
30 : !
31 : !.....-----------------------------------------------------------------
32 : USE m_constants
33 : USE m_types
34 : IMPLICIT NONE
35 : TYPE(t_cell),INTENT(IN) :: cell
36 : ! ..
37 : ! .. Scalar arguments ..
38 : INTEGER, INTENT (IN) :: xmax1,xmax2,xmax3
39 : ! ..
40 : ! .. Array arguments ..
41 : REAL, INTENT(IN) :: ro(0:xmax1*xmax2*xmax3-1)
42 : ! ..
43 : ! .. Array output ..
44 : REAL, INTENT(OUT) :: dro(0:xmax1*xmax2*xmax3-1,3)
45 : ! ..
46 : ! .. Locals ..
47 : INTEGER :: xmax(3)
48 : INTEGER :: direction,xyz(3),x1,x2,x3,i,ii(-3:2)
49 992 : REAL :: dx,drointern(0:xmax1*xmax2*xmax3-1,3)
50 :
51 :
52 496 : xmax(1)= xmax1
53 496 : xmax(2)= xmax2
54 496 : xmax(3)= xmax3
55 1984 : DO i=1,3
56 496 : IF ( xmax(i) < 3 ) THEN
57 0 : CALL juDFT_error("grid to small",calledby="grdrsis")
58 : END IF
59 : END DO
60 : IF ( (ndvgrd < 3) .or. (ndvgrd > 6) ) THEN
61 : CALL juDFT_error("ndvgrd notin [3,6]",calledby="grdrsis")
62 : ENDIF
63 :
64 :
65 3472 : DO direction=1,3
66 :
67 1488 : dx= 1./REAL(xmax(direction))
68 :
69 25792 : DO x1=0,xmax(1)-1
70 787152 : DO x2=0,xmax(2)-1
71 15150336 : DO x3=0,xmax(3)-1
72 :
73 95846400 : DO i= -3,2
74 44236800 : xyz(1)= x1
75 44236800 : xyz(2)= x2
76 44236800 : xyz(3)= x3
77 44236800 : xyz(direction)= xyz(direction)+i
78 : ! make use of periodic boundary cond. in interstitial:
79 44236800 : IF ( xyz(direction) < 0 ) THEN
80 2605056 : xyz(direction)= xyz(direction)+xmax(direction)
81 : END IF
82 44236800 : IF ( xyz(direction) >= xmax(direction) ) THEN
83 1302528 : xyz(direction)= xyz(direction)-xmax(direction)
84 : END IF
85 : ! find coordinates in 1-dim array ro:
86 51609600 : ii(i)= xyz(3)*xmax(1)*xmax(2) + xyz(2)*xmax(1) + xyz(1)
87 : END DO
88 :
89 380928 : IF (ndvgrd.EQ.3) THEN
90 : drointern(ii(0),direction)= &
91 : & df3( ro(ii(-1)), &
92 : & ro(ii(0)),ro(ii(1)), dx)
93 : ELSEIF (ndvgrd.EQ.4) THEN
94 : drointern(ii(0),direction)= &
95 : & df4( ro(ii(-1)),&
96 : & ro(ii(0)),ro(ii(1)),ro(ii(2)), dx)
97 : ELSEIF (ndvgrd.EQ.5) THEN
98 : drointern(ii(0),direction)= &
99 : & df5( ro(ii(-2)),ro(ii(-1)),&
100 : & ro(ii(0)),ro(ii(1)),ro(ii(2)), dx)
101 51609600 : ELSEIF (ndvgrd.EQ.6) THEN
102 : drointern(ii(0),direction)= &
103 : & df6( ro(ii(-3)),ro(ii(-2)),ro(ii(-1)),&
104 51609600 : & ro(ii(0)),ro(ii(1)),ro(ii(2)), dx)
105 : ENDIF
106 :
107 : END DO
108 : END DO
109 : END DO
110 :
111 : END DO
112 :
113 :
114 4915696 : DO i=0,xmax(1)*xmax(2)*xmax(3)-1
115 :
116 17203696 : DO direction=1,3
117 : dro(i,direction)= cell%bmat(1,direction)*drointern(i,1) &
118 : & + cell%bmat(2,direction)*drointern(i,2)&
119 7372800 : & + cell%bmat(3,direction)*drointern(i,3)
120 9830400 : dro(i,direction)= dro(i,direction)/(2.*pi_const)
121 : END DO
122 :
123 : END DO
124 :
125 496 : END SUBROUTINE grdrsis
126 : !--------------------------------------------------------------------
127 : ! Functions: formulae for 1st deriv.:
128 : !
129 : REAL FUNCTION df3(g1,f0,f1,d) ! three point formula
130 : REAL g1,f0,f1,d
131 : df3 = (-1*g1-0*f0+f1)/ (2*d)
132 : END FUNCTION df3
133 :
134 : REAL FUNCTION df4(g1,f0,f1,f2,d) ! four point formula
135 : REAL g1,f0,f1,f2,d
136 : df4 = (-2*g1-3*f0+6*f1-f2)/ (6*d)
137 : END FUNCTION df4
138 :
139 : REAL FUNCTION df5(g2,g1,f0,f1,f2,d) ! five point formula
140 : REAL g2,g1,f0,f1,f2,d
141 : df5 = (2*g2-16*g1-0*f0+16*f1-2*f2)/ (24*d)
142 : END FUNCTION df5
143 :
144 : REAL FUNCTION df6(g3,g2,g1,f0,f1,f2,d) ! six point formula
145 : REAL g3,g2,g1,f0,f1,f2,d
146 7372800 : df6 = (-4*g3+30*g2-120*g1+40*f0+60*f1-6*f2)/ (120*d)
147 : END FUNCTION df6
148 :
149 : !----------------------------------------------------------------------
150 : END MODULE m_grdrsis
|
