|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00010 MODULE ps_wavelet_scaling_function 00011 USE kinds, ONLY: dp 00012 USE lazy, ONLY: lazy_arrays 00013 #include "cp_common_uses.h" 00014 00015 PRIVATE 00016 00017 PUBLIC :: scaling_function,& 00018 scf_recursion 00019 00020 CONTAINS 00021 00022 ! ***************************************************************************** 00025 SUBROUTINE scaling_function(itype,nd,nrange,a,x) 00026 00027 !Type of interpolating functions 00028 INTEGER, INTENT(in) :: itype, nd 00029 INTEGER, INTENT(out) :: nrange 00030 REAL(KIND=dp), DIMENSION(0:nd), 00031 INTENT(out) :: a, x 00032 00033 INTEGER :: i, i_all, m, ni, nt 00034 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: y 00035 REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht 00036 00037 !Number of points: must be 2**nex 00038 00039 a = 0.0_dp 00040 x = 0.0_dp 00041 m=itype+2 00042 CALL lazy_arrays(itype,m,ch,cg,cgt,cht) 00043 00044 ni=2*itype 00045 nrange = ni 00046 ALLOCATE(y(0:nd),stat=i_all) 00047 IF (i_all /= 0) THEN 00048 WRITE(*,*)' scaling_function: problem of memory allocation' 00049 STOP 00050 END IF 00051 00052 ! plot scaling function 00053 CALL zero(nd+1,x) 00054 CALL zero(nd+1,y) 00055 nt=ni 00056 x(nt/2-1)=1._dp 00057 loop1: DO 00058 nt=2*nt 00059 00060 CALL back_trans(nd,nt,x,y,m,ch,cg,cgt,cht) 00061 CALL dcopy(nt,y,1,x,1) 00062 IF (nt.EQ.nd) THEN 00063 EXIT loop1 00064 END IF 00065 END DO loop1 00066 00067 !open (unit=1,file='scfunction',status='unknown') 00068 DO i=0,nd 00069 a(i) = 1._dp*i*ni/nd-(.5_dp*ni-1._dp) 00070 !write(1,*) 1._dp*i*ni/nd-(.5_dp*ni-1._dp),x(i) 00071 END DO 00072 !close(1) 00073 DEALLOCATE(ch,cg,cgt,cht) 00074 DEALLOCATE(y,stat=i_all) 00075 IF (i_all /= 0) THEN 00076 WRITE(*,*)' scaling_function: problem of memory deallocation' 00077 STOP 00078 END IF 00079 END SUBROUTINE scaling_function 00080 00081 ! ***************************************************************************** 00084 SUBROUTINE wavelet_function(itype,nd,a,x) 00085 00086 !Type of the interpolating scaling function 00087 INTEGER, INTENT(in) :: itype, nd 00088 REAL(KIND=dp), DIMENSION(0:nd), 00089 INTENT(out) :: a, x 00090 00091 INTEGER :: i, i_all, m, ni, nt 00092 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: y 00093 REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht 00094 00095 !must be 2**nex 00096 00097 a = 0.0_dp 00098 x = 0.0_dp 00099 m=itype+2 00100 ni=2*itype 00101 CALL lazy_arrays(itype,m,ch,cg,cgt,cht) 00102 ALLOCATE(y(0:nd),stat=i_all) 00103 IF (i_all /= 0) THEN 00104 WRITE(*,*)' wavelet_function: problem of memory allocation' 00105 STOP 00106 END IF 00107 00108 ! plot wavelet 00109 CALL zero(nd+1,x) 00110 CALL zero(nd+1,y) 00111 nt=ni 00112 x(nt+nt/2-1)=1._dp 00113 loop3: DO 00114 nt=2*nt 00115 !write(6,*) 'nd,nt',nd,nt 00116 CALL back_trans(nd,nt,x,y,m,ch,cg,cgt,cht) 00117 CALL dcopy(nd,y,1,x,1) 00118 IF (nt.EQ.nd) THEN 00119 EXIT loop3 00120 END IF 00121 END DO loop3 00122 00123 !open (unit=1,file='wavelet',status='unknown') 00124 DO i=0,nd-1 00125 a(i) = 1._dp*i*ni/nd-(.5_dp*ni-.5_dp) 00126 !write(1,*) 1._dp*i*ni/nd-(.5_dp*ni-.5_dp),x(i) 00127 END DO 00128 !close(1) 00129 DEALLOCATE(ch,cg,cgt,cht) 00130 DEALLOCATE(y,stat=i_all) 00131 IF (i_all /= 0) THEN 00132 WRITE(*,*)' wavelet_function: problem of memory deallocation' 00133 STOP 00134 END IF 00135 00136 END SUBROUTINE wavelet_function 00137 00138 ! ***************************************************************************** 00142 SUBROUTINE scf_recursion(itype,n_iter,n_range,kernel_scf,kern_1_scf) 00143 INTEGER, INTENT(in) :: itype, n_iter, n_range 00144 REAL(KIND=dp), INTENT(inout) :: kernel_scf(-n_range:n_range) 00145 REAL(KIND=dp), INTENT(out) :: kern_1_scf(-n_range:n_range) 00146 00147 INTEGER :: m 00148 REAL(KIND=dp), DIMENSION(:), POINTER :: cg, cgt, ch, cht 00149 00150 kern_1_scf = 0.0_dp 00151 m=itype+2 00152 CALL lazy_arrays(itype,m,ch,cg,cgt,cht) 00153 CALL scf_recurs(n_iter,n_range,kernel_scf,kern_1_scf,m,ch,cg,cgt,cht) 00154 DEALLOCATE(ch,cg,cgt,cht) 00155 00156 END SUBROUTINE scf_recursion 00157 00158 ! ***************************************************************************** 00161 SUBROUTINE zero(n,x) 00162 INTEGER, INTENT(in) :: n 00163 REAL(KIND=dp), INTENT(out) :: x(n) 00164 00165 INTEGER :: i 00166 00167 DO i=1,n 00168 x(i)=0._dp 00169 END DO 00170 END SUBROUTINE zero 00171 00172 ! ***************************************************************************** 00179 SUBROUTINE for_trans(nd,nt,x,y,m,ch,cg,cgt,cht) 00180 INTEGER, INTENT(in) :: nd, nt 00181 REAL(KIND=dp), INTENT(in) :: x(0:nd-1) 00182 REAL(KIND=dp), INTENT(out) :: y(0:nd-1) 00183 INTEGER :: m 00184 REAL(KIND=dp), DIMENSION(:), POINTER :: ch, cg, cgt, cht 00185 00186 INTEGER :: i, ind, j 00187 00188 y = 0.0_dp 00189 DO i=0,nt/2-1 00190 y( i)=0._dp 00191 y(nt/2+i)=0._dp 00192 00193 DO j=-m+1,m 00194 00195 ! periodically wrap index if necessary 00196 ind=j+2*i 00197 loop99: DO 00198 IF (ind.lt.0) THEN 00199 ind=ind+nt 00200 CYCLE loop99 00201 END IF 00202 IF (ind.ge.nt) THEN 00203 ind=ind-nt 00204 CYCLE loop99 00205 END IF 00206 EXIT loop99 00207 END DO loop99 00208 00209 y( i)=y( i)+cht(j)*x(ind) 00210 y(nt/2+i)=y(nt/2+i)+cgt(j)*x(ind) 00211 END DO 00212 00213 END DO 00214 00215 END SUBROUTINE for_trans 00216 00217 ! ***************************************************************************** 00218 SUBROUTINE back_trans(nd,nt,x,y,m,ch,cg,cgt,cht) 00219 ! backward wavelet transform 00220 ! nd: length of data set 00221 ! nt length of data in data set to be transformed 00222 ! m filter length (m has to be even!) 00223 ! x input data, y output data 00224 INTEGER, INTENT(in) :: nd, nt 00225 REAL(KIND=dp), INTENT(in) :: x(0:nd-1) 00226 REAL(KIND=dp), INTENT(out) :: y(0:nd-1) 00227 INTEGER :: m 00228 REAL(KIND=dp), DIMENSION(:), POINTER :: ch, cg, cgt, cht 00229 00230 INTEGER :: i, ind, j 00231 00232 y = 0.0_dp 00233 00234 DO i=0,nt/2-1 00235 y(2*i+0)=0._dp 00236 y(2*i+1)=0._dp 00237 00238 DO j=-m/2,m/2-1 00239 00240 ! periodically wrap index if necessary 00241 ind=i-j 00242 loop99: DO 00243 IF (ind.lt.0) THEN 00244 ind=ind+nt/2 00245 CYCLE loop99 00246 END IF 00247 IF (ind.ge.nt/2) THEN 00248 ind=ind-nt/2 00249 CYCLE loop99 00250 END IF 00251 EXIT loop99 00252 END DO loop99 00253 00254 y(2*i+0)=y(2*i+0) + ch(2*j-0)*x(ind)+cg(2*j-0)*x(ind+nt/2) 00255 y(2*i+1)=y(2*i+1) + ch(2*j+1)*x(ind)+cg(2*j+1)*x(ind+nt/2) 00256 END DO 00257 00258 END DO 00259 00260 END SUBROUTINE back_trans 00261 00262 ! ***************************************************************************** 00265 SUBROUTINE ftest(m,ch,cg,cgt,cht) 00266 INTEGER :: m 00267 REAL(KIND=dp), DIMENSION(:), POINTER :: ch, cg, cgt, cht 00268 00269 CHARACTER(len=*), PARAMETER :: fmt22 = "(a,i3,i4,4(e17.10))" 00270 00271 INTEGER :: i, j, l 00272 REAL(KIND=dp) :: eps, t1, t2, t3, t4 00273 00274 ! do i=-m,m 00275 ! write(6,*) i,ch(i),cg(i) 00276 ! end do 00277 00278 DO i=-m,m 00279 DO j=-m,m 00280 t1=0._dp 00281 t2=0._dp 00282 t3=0._dp 00283 t4=0._dp 00284 DO l=-3*m,3*m 00285 IF ( l-2*i.ge.-m .AND. l-2*i.le.m .AND. & 00286 l-2*j.ge.-m .AND. l-2*j.le.m ) THEN 00287 t1=t1+ch(l-2*i)*cht(l-2*j) 00288 t2=t2+cg(l-2*i)*cgt(l-2*j) 00289 t3=t3+ch(l-2*i)*cgt(l-2*j) 00290 t4=t4+cht(l-2*i)*cg(l-2*j) 00291 END IF 00292 END DO 00293 eps=1.e-10_dp 00294 IF (i.eq.j) THEN 00295 IF (ABS(t1-1._dp).gt.eps .OR. ABS(t2-1._dp).gt.eps .OR. & 00296 & ABS(t3).gt.eps .OR. ABS(t4).gt.eps ) THEN 00297 WRITE(6,fmt22) 'Orthogonality ERROR', i,j,t1,t2,t3,t4 00298 END IF 00299 ELSE 00300 IF (ABS(t1).gt.eps .OR. ABS(t2).gt.eps .OR. & 00301 & ABS(t3).gt.eps .OR. ABS(t4).gt.eps ) THEN 00302 WRITE(6,fmt22) 'Orthogonality ERROR', i,j,t1,t2,t3,t4 00303 END IF 00304 END IF 00305 END DO 00306 END DO 00307 00308 WRITE(6,*) 'FILTER TEST PASSED' 00309 00310 END SUBROUTINE ftest 00311 00312 ! ***************************************************************************** 00316 SUBROUTINE scf_recurs(n_iter,n_range,kernel_scf,kern_1_scf,m,ch,cg,cgt,cht) 00317 INTEGER, INTENT(in) :: n_iter, n_range 00318 REAL(KIND=dp), INTENT(inout) :: kernel_scf(-n_range:n_range) 00319 REAL(KIND=dp), INTENT(out) :: kern_1_scf(-n_range:n_range) 00320 INTEGER :: m 00321 REAL(KIND=dp), DIMENSION(:), POINTER :: ch, cg, cgt, cht 00322 00323 INTEGER :: i, i_iter, ind, j 00324 REAL(KIND=dp) :: kern, kern_tot 00325 00326 kern_1_scf = 0.0_dp 00327 !Start the iteration to go from p0gauss to pgauss 00328 loop_iter_scf: DO i_iter=1,n_iter 00329 kern_1_scf(:) = kernel_scf(:) 00330 kernel_scf(:) = 0._dp 00331 loop_iter_i: DO i=0,n_range 00332 kern_tot = 0._dp 00333 DO j=-m,m 00334 ind = 2*i-j 00335 IF (ABS(ind) > n_range) THEN 00336 kern = 0._dp 00337 ELSE 00338 kern = kern_1_scf(ind) 00339 END IF 00340 kern_tot = kern_tot + ch(j)*kern 00341 END DO 00342 IF (kern_tot == 0._dp) THEN 00343 !zero after (be sure because strictly == 0._dp) 00344 EXIT loop_iter_i 00345 ELSE 00346 kernel_scf( i) = 0.5_dp*kern_tot 00347 kernel_scf(-i) = kernel_scf(i) 00348 END IF 00349 END DO loop_iter_i 00350 END DO loop_iter_scf 00351 END SUBROUTINE scf_recurs 00352 00353 END MODULE ps_wavelet_scaling_function
1.7.3