TABLE OF CONTENTS


ABINIT/psp3nl [ Functions ]

[ Top ] [ Functions ]

NAME

 psp3nl

FUNCTION

 Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998).
 Uses Gaussians for fully nonlocal form, analytic expressions.

COPYRIGHT

 Copyright (C) 1998-2017 ABINIT group (DCA, FRD, XG, GMR, PT)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  h11s=factor defining strength of 1st projector for l=0 channel
  h22s=factor defining strength of 2nd projector for l=0 channel
  h33s=factor defining strength of 3rd projector for l=0 channel
  h11p=factor defining strength of 1st projector for l=1 channel
  h22p=factor defining strength of 2nd projector for l=1 channel
  h33p=factor defining strength of 2nd projector for l=1 channel
  h11d=factor defining strength of 1st projector for l=2 channel
  h22d=factor defining strength of 2nd projector for l=2 channel
  h33d=factor defining strength of 2nd projector for l=2 channel
  h11f=factor defining strength of 1st projector for l=3 channel
  mproj=maximum number of projectors in any channel
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mqgrid=number of grid points for qgrid
  qgrid(mqgrid)=array of |G| values
  rrd=core radius for d channel (bohr)
  rrf=core radius for f channel (bohr)
  rrp=core radius for p channel (bohr)
  rrs=core radius for s channel (bohr)

OUTPUT

  ekb(mpsang,mproj)=Kleinman-Bylander energies
  ffspl(mqgrid,2,mpssang,mproj)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum and
   each projectors

PARENTS

      psp3in

CHILDREN

      spline,zhpev

SOURCE

 51 #if defined HAVE_CONFIG_H
 52 #include "config.h"
 53 #endif
 54 
 55 #include "abi_common.h"
 56 
 57 subroutine psp3nl(ekb,ffspl,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,&
 58 &                  h33d,h11f,mproj,mpsang,mqgrid,qgrid,rrd,rrf,rrp,rrs)
 59 
 60  use defs_basis
 61  use m_splines
 62  use m_errors
 63  use m_profiling_abi
 64 
 65 !This section has been created automatically by the script Abilint (TD).
 66 !Do not modify the following lines by hand.
 67 #undef ABI_FUNC
 68 #define ABI_FUNC 'psp3nl'
 69 !End of the abilint section
 70 
 71  implicit none
 72 
 73 !Arguments ------------------------------------
 74 !scalars
 75  integer,intent(in) :: mproj,mpsang,mqgrid
 76  real(dp),intent(in) :: h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s,rrd
 77  real(dp),intent(in) :: rrf,rrp,rrs
 78 !arrays
 79  real(dp),intent(in) :: qgrid(mqgrid)
 80  real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj)
 81 
 82 !Local variables-------------------------------
 83 !scalars
 84  integer :: info,iproj,iqgrid,ldz,mu,nproj,nu
 85  real(dp) :: qmax
 86  character(len=500) :: message
 87  character :: jobz,uplo
 88 !arrays
 89  real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3)
 90  real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:)
 91 
 92 ! *************************************************************************
 93 
 94 !DEBUG
 95 !write(std_out,*)' psp3nl : enter '
 96 !stop
 97 !ENDDEBUG
 98 
 99  ABI_ALLOCATE(ppspl,(mqgrid,2,mpsang,mproj))
100  ABI_ALLOCATE(work,(mqgrid))
101 
102  qmax=qgrid(mqgrid)
103  jobz='v'
104  uplo='u'
105  ekb(:,:)=0.0d0
106 
107 !---------------------------------------------------------------
108 !Treat s channel
109 
110  nproj=0
111  ap(:,:)=0.0d0
112 !If there is at least one s-projector
113  if  ( abs(h11s) >= 1.0d-8 ) then
114    nproj=1 ; ldz=1 ; ap(1,1)=h11s
115  end if
116  nproj=1
117 !If there is a second projector
118  if  ( abs(h22s) >= 1.0d-8 ) then
119    nproj=2 ; ldz=2 ; ap(1,3)=h22s
120    ap(1,2)=-0.5d0*sqrt(0.6d0)*h22s
121  end if
122 !If there is a third projector
123  if ( abs(h33s) >= 1.0d-8 ) then
124    nproj=3 ; ldz=3 ; ap(1,6)=h33s
125    ap(1,4)=0.5d0*sqrt(5.d0/21.d0)*h33s
126    ap(1,5)=-0.5d0*sqrt(100.d0/63.d0)*h33s
127  end if
128 
129  if(nproj/=0)then
130 
131    ABI_ALLOCATE(uu,(nproj,nproj))
132    ABI_ALLOCATE(zz,(2,nproj,nproj))
133 
134    if (nproj > 1) then
135      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
136      uu(:,:)=zz(1,:,:)
137    else
138      ww(1)=h11s
139      uu(1,1)=1.0d0
140    end if
141 
142 !  Initialization of ekb, and spline fitting
143    do iproj=1,nproj
144      ekb(1,iproj)=ww(iproj)*32.d0*(rrs**3)*(pi**2.5d0)/(4.d0*pi)**2
145      if(iproj==1)then
146        do iqgrid=1,mqgrid
147          ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2)
148        end do
149        yp1j(1)=0.d0
150        ypnj(1)=-(two_pi*rrs)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrs)**2)
151      else if(iproj==2)then
152        do iqgrid=1,mqgrid
153          ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0)     &
154 &         *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) &
155 &         *( 3.d0-(two_pi*qgrid(iqgrid)*rrs)**2 )
156        end do
157        yp1j(2)=0.0d0
158        ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrs)**2*qmax &
159 &       *exp(-0.5d0*(two_pi*qmax*rrs)**2) * (-5.d0+(two_pi*qmax*rrs)**2)
160      else if(iproj==3)then
161        do iqgrid=1,mqgrid
162          ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*&
163 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * &
164 &         (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrs)**2 + &
165 &         (two_pi*qgrid(iqgrid)*rrs)**4)
166        end do
167        yp1j(3)=0.0d0
168        ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrs)**2) * &
169 &       (two_pi*rrs)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrs)**2-(two_pi*qmax*rrs)**4)
170      end if
171      call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,&
172 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj))
173    end do
174 
175 !  Linear combination using the eigenvectors
176    ffspl(:,:,1,:)=0.0d0
177    do mu=1,nproj
178      do nu=1,nproj
179        do iqgrid=1,mqgrid
180          ffspl(iqgrid,1:2,1,mu)=ffspl(iqgrid,1:2,1,mu) &
181 &         +uu(nu,mu)*ppspl(iqgrid,1:2,1,nu)
182        end do
183      end do
184    end do
185 
186    ABI_DEALLOCATE(uu)
187    ABI_DEALLOCATE(zz)
188 
189 !  End condition on nproj(/=0)
190  end if
191 
192 !DEBUG
193 !write(std_out,*)' psp3nl : after s channel '
194 !stop
195 !ENDDEBUG
196 
197 !--------------------------------------------------------------------
198 !Now treat p channel
199 
200  nproj=0
201  ap(:,:)=0.0d0
202 !If there is at least one projector
203  if  ( abs(h11p) >= 1.0d-8 ) then
204    nproj=1 ; ldz=1 ; ap(1,1)=h11p
205  end if
206 !If there is a second projector
207  if  ( abs(h22p) >= 1.0d-8 ) then
208    nproj=2 ; ldz=2 ; ap(1,3)=h22p
209    ap(1,2)=-0.5d0*sqrt(5.d0/7.d0)*h22p
210  end if
211 !If there is a third projector
212  if ( abs(h33p) >= 1.0d-8 ) then
213    nproj=3 ; ldz=3 ; ap(1,6)=h33p
214    ap(1,4)= (1.d0/6.d0)*sqrt(35.d0/11.d0)*h33p
215    ap(1,5)=-(1.d0/6.d0)*(14.d0/sqrt(11.d0))*h33p
216  end if
217 
218  if(nproj/=0)then
219 
220    ABI_ALLOCATE(uu,(nproj,nproj))
221    ABI_ALLOCATE(zz,(2,nproj,nproj))
222 
223    if (nproj > 1) then
224      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
225      uu(:,:)=zz(1,:,:)
226    else
227      ww(1)=h11p
228      uu(1,1)=1.0d0
229    end if
230 
231 !  Initialization of ekb, and spline fitting
232    do iproj=1,nproj
233      ekb(2,iproj)=ww(iproj)*64.d0*(rrp**5)*(pi**2.5d0)/(4.d0*pi)**2
234      if(iproj==1)then
235        do iqgrid=1,mqgrid
236          ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* &
237 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * (two_pi*qgrid(iqgrid))
238        end do
239        yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0))
240        ypnj(1)=-two_pi*((two_pi*qmax*rrp)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrp)**2)*&
241 &       (1.0d0/sqrt(3.0d0))
242      else if(iproj==2)then
243        do iqgrid=1,mqgrid
244          ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* &
245 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
246 &         (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrp)**2)
247        end do
248        yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0))
249        ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
250 &       (-8*(two_pi*qmax*rrp)**2 + (two_pi*qmax*rrp)**4 + 5.0d0)
251      else if(iproj==3)then
252        do iqgrid=1,mqgrid
253          ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*&
254 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
255 &         (two_pi*qgrid(iqgrid))*&
256 &         (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrp)**2+(two_pi*qgrid(iqgrid)*rrp)**4)
257        end do
258        yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0)
259        ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
260 &       (35.0d0-77.0d0*(two_pi*qmax*rrp)**2+19.0d0*(two_pi*qmax*rrp)**4 - &
261 &       (two_pi*qmax*rrp)**6)
262      end if
263      call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,&
264 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj))
265    end do
266 
267 !  Linear combination using the eigenvectors
268    ffspl(:,:,2,:)=0.0d0
269    do mu=1,nproj
270      do nu=1,nproj
271        do iqgrid=1,mqgrid
272          ffspl(iqgrid,1:2,2,mu)=ffspl(iqgrid,1:2,2,mu) &
273 &         +uu(nu,mu)*ppspl(iqgrid,1:2,2,nu)
274        end do
275      end do
276    end do
277 
278    ABI_DEALLOCATE(uu)
279    ABI_DEALLOCATE(zz)
280 
281 !  End condition on nproj(/=0)
282  end if
283 
284 !DEBUG
285 !write(std_out,*)' psp3nl : after p channel '
286 !stop
287 !ENDDEBUG
288 
289 !-----------------------------------------------------------------------
290 !Now treat d channel.
291 
292  nproj=0
293  ap(:,:)=0.0d0
294 !If there is at least one projector
295  if  ( abs(h11d) >= 1.0d-8 ) then
296    nproj=1 ; ldz=1 ; ap(1,1)=h11d
297  end if
298 !If there is a second projector
299  if  ( abs(h22d) >= 1.0d-8 ) then
300    nproj=2 ; ldz=2 ; ap(1,3)=h22d
301    ap(1,2)=-0.5d0*sqrt(7.d0/9.d0)*h22d
302  end if
303 !If there is a third projector. Warning : only two projectors are allowed.
304  if ( abs(h33d) >= 1.0d-8 ) then
305    write(message, '(a,a,a)' )&
306 &   '  only two d-projectors are allowed ',ch10,&
307 &   '  Action : check your pseudopotential file.'
308    MSG_ERROR(message)
309 !  nproj=3 ; ldz=3 ; ap(1,6)=h33d
310 !  ap(1,4)= 0.5d0*sqrt(63.d0/143.d0)*h33d
311 !  ap(1,5)= -0.5d0*(18.d0/sqrt(143.d0))*h33d
312  end if
313 
314  if(nproj/=0)then
315 
316    ABI_ALLOCATE(uu,(nproj,nproj))
317    ABI_ALLOCATE(zz,(2,nproj,nproj))
318 
319    if (nproj > 1) then
320      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
321      uu(:,:)=zz(1,:,:)
322    else
323      ww(1)=h11d
324      uu(1,1)=1.0d0
325    end if
326 
327 !  Initialization of ekb, and spline fitting
328    do iproj=1,nproj
329      ekb(3,iproj)=ww(iproj)*128.d0*(rrd**7)*(pi**2.5d0)/(4.d0*pi)**2
330      if(iproj==1)then
331        do iqgrid=1,mqgrid
332          ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* &
333 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * (two_pi*qgrid(iqgrid))**2
334        end do
335        yp1j(1)=0.0d0
336        ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*&
337 &       exp(-0.5d0*(two_pi*qmax*rrd)**2)*qmax*(2d0-(two_pi*qmax*rrd)**2)
338      else if(iproj==2)then
339        do iqgrid=1,mqgrid
340          ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* &
341 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * &
342 &         ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrd)**2)
343        end do
344        yp1j(2)=0.0d0
345        ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrd)**2)* &
346 &       qmax*(two_pi**2)*( (two_pi*qmax*rrd)**4 - 11.0d0*(two_pi*qmax*rrd)**2 + 14.0d0)
347      end if
348      call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,&
349 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj))
350    end do
351 
352 !  Linear combination using the eigenvectors
353    ffspl(:,:,3,:)=0.0d0
354    do mu=1,nproj
355      do nu=1,nproj
356        do iqgrid=1,mqgrid
357          ffspl(iqgrid,1:2,3,mu)=ffspl(iqgrid,1:2,3,mu) &
358 &         +uu(nu,mu)*ppspl(iqgrid,1:2,3,nu)
359        end do
360      end do
361    end do
362 
363    ABI_DEALLOCATE(uu)
364    ABI_DEALLOCATE(zz)
365 
366 !  End condition on nproj(/=0)
367  end if
368 
369 !DEBUG
370 !write(std_out,*)' psp3nl : after d channel '
371 !stop
372 !ENDDEBUG
373 
374 !-----------------------------------------------------------------------
375 !Treat now f channel (max one projector ! - so do not use ppspl)
376 
377 !l=3 first projector
378  if (abs(h11f)>1.d-12) then
379    ekb(4,1)=h11f*(256.0d0/105.0d0)*(rrf**9)*(pi**2.5d0)/(4.d0*pi)**2
380    do iqgrid=1,mqgrid
381      ffspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* &
382 &     exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrf)**2)
383    end do
384 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
385    yp1j(1)=0d0
386    ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrf)**2)*&
387 &   (3.0d0-(two_pi*qmax*rrf)**2)
388 !  Fit spline to get second derivatives by spline fit
389    call spline(qgrid,ffspl(:,1,4,1),mqgrid,&
390 &   yp1j(1),ypnj(1),ffspl(:,2,4,1))
391  end if
392 
393 !-----------------------------------------------------------------------
394 
395  ABI_DEALLOCATE(ppspl)
396  ABI_DEALLOCATE(work)
397 
398 !DEBUG
399 !write(std_out,*)' psp3nl : exit '
400 !stop
401 !ENDDEBUG
402 
403 end subroutine psp3nl