TABLE OF CONTENTS


ABINIT/getgh2c [ Functions ]

[ Top ] [ Functions ]

NAME

 getgh2c

FUNCTION

 Compute <G|H^(2)|C> (or <G|H^(2)-Eps.S^(2)|C>) for input vector |C> expressed in reciprocal space.
 (H^(2) is the 2nd-order pertubed Hamiltonian, S^(2) is the 2nd-order perturbed overlap operator).
 Result is put in array gh2c.
 If required, part of <G|K(2)+Vnonlocal^(2)|C> not depending on VHxc^(2) is also returned in gvnl2.
 If required, <G|S^(2)|C> is returned in gs2c (S=overlap - PAW only)
 Available for the following cases :
  ipert = natom+10 (dkdk)   :       2nd derivative w.r.t wavevector
          natom+11 (dkdE)   : mixed 2nd derivative w.r.t wavector     and eletric field
  also if natom+12<=ipert<=2*natom+11 :
                   (dtaudE) : mixed 2nd derivative w.r.t atom. displ. and eletric field (nonlocal only)

INPUTS

  cwavef(2,npw*nspinor)=input wavefunction, in reciprocal space
  cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C>
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian
  idir=direction of the perturbation
  ipert=type of the perturbation
  lambda=real use to apply H^(2)-lambda.S^(2)
  mpi_enreg=information about MPI parallelization
  optlocal=0: local part of H^(2) is not computed in gh2c=<G|H^(2)|C>
           1: local part of H^(2) is computed in gh2c=<G|H^(2)|C>
  optnl=0: non-local part of H^(2) is not computed in gh2c=<G|H^(2)|C>
        1: non-local part of H^(2) depending on VHxc^(2) is not computed in gh2c=<G|H^(2)|C>
        2: non-local part of H^(2) is totally computed in gh2c=<G|H^(2)|C>
  opt_gvnl2=option controlling the use of gvnl2 array:
            0: not used
            1: used as input:    - used only for PAW and ipert=natom+11/+12
               At input: contains the derivative w.r.t wavevector of cwavef (times i)
  rf_hamkq <type(rf_hamiltonian_type)>=all data for the 2nd-order Hamiltonian at k,k+q
  sij_opt= -PAW ONLY-  if  0, only matrix elements <G|H^(2)|C> have to be computed
     (S=overlap)       if  1, matrix elements <G|S^(2)|C> have to be computed in gs2c in addition to gh2c
                       if -1, matrix elements <G|H^(2)-lambda.S^(2)|C> have to be computed in gh2c (gs2c not used)
  tim_getgh2c=timing code of the calling subroutine (can be set to 0 if not attributed)
  usevnl=1 if gvnl2=(part of <G|K^(2)+Vnl^(2)-lambda.S^(2)|C> not depending on VHxc^(2)) has to be input/output

OUTPUT

 gh2c(2,npw1*nspinor)= <G|H^(2)|C> or  <G|H^(2)-lambda.S^(2)|C>
                     (only kinetic+non-local parts if optlocal=0)
 if (usevnl==1)
  gvnl2(2,npw1*nspinor)=  part of <G|K^(2)+Vnl^(2)|C> not depending on VHxc^(2)              (sij_opt/=-1)
                       or part of <G|K^(2)+Vnl^(2)-lambda.S^(2)|C> not depending on VHxc^(2) (sij_opt==-1)
 if (sij_opt=1)
  gs2c(2,npw1*nspinor)=<G|S^(2)|C> (S=overlap).

PARENTS

      m_rf2

CHILDREN

      nonlop,pawcprj_alloc,pawcprj_free

SOURCE

106 subroutine getgh2c(cwavef,cwaveprj,gh2c,gs2c,gs_hamkq,gvnl2,idir,ipert,lambda,&
107 &                  mpi_enreg,optlocal,optnl,opt_gvnl2,rf_hamkq,sij_opt,tim_getgh2c,usevnl,conj,enl,optkin)
108 
109 
110 !This section has been created automatically by the script Abilint (TD).
111 !Do not modify the following lines by hand.
112 #undef ABI_FUNC
113 #define ABI_FUNC 'getgh2c'
114 !End of the abilint section
115 
116  implicit none
117 
118 !Arguments ------------------------------------
119 !scalars
120  logical,intent(in),optional :: conj
121  integer,intent(in) :: idir,ipert,optlocal,optnl,opt_gvnl2,sij_opt,tim_getgh2c,usevnl
122  integer,intent(in),optional :: optkin
123  real(dp),intent(in) :: lambda
124  type(MPI_type),intent(in) :: mpi_enreg
125  type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq
126  type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq
127 !arrays
128  real(dp),intent(in),optional,target :: enl(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,gs_hamkq%dimekbq)
129  real(dp),intent(inout) :: cwavef(:,:)
130  real(dp),intent(inout),target :: gvnl2(:,:)
131  real(dp),intent(out) :: gh2c(:,:),gs2c(:,:)
132  type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
133 
134 !Local variables-------------------------------
135 !scalars
136  integer,parameter :: tim_nonlop=0
137  integer :: choice,cpopt,iatm,idir1,idir2,idirc,ipw,ipws,ispinor,my_nspinor
138  integer :: natom,ncpgr,nnlout=1,npw,npw1,paw_opt,signs,usecprj
139  logical :: compute_conjugate,has_kin,has_vnl,pert_phon_elfd
140  real(dp) :: enlout_dum(1)
141  character(len=500) :: msg
142 !arrays
143 ! integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
144  integer,parameter :: alpha(9)=(/1,2,3,2,1,1,3,3,2/),beta(9)=(/1,2,3,3,3,2,2,1,1/)
145 !real(dp) :: tsec(2)
146  real(dp) :: svectout_dum(1,1),vectout_dum(1,1)
147  real(dp),allocatable :: nonlop_out(:,:)
148  real(dp), pointer :: gvnl2_(:,:)
149  real(dp), pointer :: ddkinpw(:),kinpw1(:),enl_ptr(:,:,:,:)
150  real(dp),allocatable,target :: enl_temp(:,:,:,:)
151  type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:)
152  type(pawcprj_type),pointer :: cwaveprj_ptr(:,:)
153 
154 ! *********************************************************************
155 
156  DBG_ENTER("COLL")
157 
158 !Keep track of total time spent in getgh2c
159 !call timab(196+tim_getgh2c,1,tsec)
160 
161 !======================================================================
162 !== Initialisations and compatibility tests
163 !======================================================================
164 
165  npw  =gs_hamkq%npw_k
166  npw1 =gs_hamkq%npw_kp
167  natom=gs_hamkq%natom
168 
169 !Compatibility tests
170  if(ipert/=natom+10.and.ipert/=natom+11.and.ipert>2*natom+11)then
171    msg='only ipert==natom+10/+11 and natom+11<=ipert<=2*natom+11 implemented!'
172    MSG_BUG(msg)
173  end if
174  pert_phon_elfd = .false.
175  if (ipert>natom+11.and.ipert<=2*natom+11) pert_phon_elfd = .true.
176  if (mpi_enreg%paral_spinor==1) then
177    msg='Not compatible with parallelization over spinorial components!'
178    MSG_BUG(msg)
179  end if
180  if (gs_hamkq%nvloc>1) then
181    msg='Not compatible with nvloc=4 (non-coll. magnetism)!'
182    MSG_BUG(msg)
183  end if
184  if((ipert==natom+11.or.pert_phon_elfd).and.gs_hamkq%usepaw==1.and.optnl>=1) then
185    if (gs_hamkq%nvloc>1) then
186      msg='Not compatible with nvloc=4 (non-coll. magnetism)!'
187      MSG_BUG(msg)
188    end if
189    if (present(enl)) then
190      enl_ptr => enl
191    else if (associated(rf_hamkq%e1kbfr).and.associated(rf_hamkq%e1kbsc).and.optnl==2) then
192      ABI_CHECK(size(rf_hamkq%e1kbfr,4)==1,'BUG in getgh2c: cplex_rf>1!')
193      ABI_CHECK(size(rf_hamkq%e1kbsc,4)==1,'BUG in getgh2c: cplex_rf>1!')
194      ABI_ALLOCATE(enl_temp,(gs_hamkq%dimekb1,gs_hamkq%dimekb2,gs_hamkq%nspinor**2,gs_hamkq%dimekbq))
195      enl_temp(:,:,:,:) = rf_hamkq%e1kbfr(:,:,:,:) + rf_hamkq%e1kbsc(:,:,:,:)
196      enl_ptr => enl_temp
197    else if (associated(rf_hamkq%e1kbfr)) then
198      ABI_CHECK(size(rf_hamkq%e1kbfr,4)==1,'BUG in getgh2c: cplex_rf>1!')
199      enl_ptr => rf_hamkq%e1kbfr
200    else
201      msg='For ipert=natom+11/pert_phon_elfd : e1kbfr and/or e1kbsc must be associated or enl optional input must be present.'
202      MSG_BUG(msg)
203    end if
204    if (usevnl==0) then
205      msg='gvnl2 must be allocated for ipert=natom+11/pert_phon_elfd !'
206      MSG_BUG(msg)
207    end if
208    if(opt_gvnl2==0) then
209      msg='opt_gvnl2=0 not compatible with ipert=natom+11/pert_phon_elfd !'
210      MSG_BUG(msg)
211    end if
212  end if
213 
214 !Check sizes
215  my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor)
216  if (size(cwavef)<2*npw*my_nspinor) then
217    msg='wrong size for cwavef!'
218    MSG_BUG(msg)
219  end if
220  if (size(gh2c)<2*npw1*my_nspinor) then
221    msg='wrong size for gh2c!'
222    MSG_BUG(msg)
223  end if
224  if (usevnl/=0) then
225    if (size(gvnl2)<2*npw1*my_nspinor) then
226      msg='wrong size for gvnl2!'
227      MSG_BUG(msg)
228    end if
229  end if
230  if (sij_opt==1) then
231    if (size(gs2c)<2*npw1*my_nspinor) then
232      msg='wrong size for gs2c!'
233      MSG_BUG(msg)
234    end if
235  end if
236 
237 !PAW: specific treatment for usecprj input arg
238 !     force it to zero if cwaveprj is not allocated
239  usecprj=gs_hamkq%usecprj ; ncpgr=0
240  if(gs_hamkq%usepaw==1) then
241    if (size(cwaveprj)==0) usecprj=0
242    if (usecprj/=0) then
243      ncpgr=cwaveprj(1,1)%ncpgr
244      if (size(cwaveprj)<gs_hamkq%natom*my_nspinor) then
245        msg='wrong size for cwaveprj!'
246        MSG_BUG(msg)
247      end if
248    end if
249  else
250    if(usecprj==1)then
251      msg='usecprj==1 not allowed for NC psps !'
252      MSG_BUG(msg)
253    end if
254  end if
255 
256 ! tim_nonlop=8
257 ! if (tim_getgh2c==1.and.ipert<=natom) tim_nonlop=7
258 ! if (tim_getgh2c==2.and.ipert<=natom) tim_nonlop=5
259 ! if (tim_getgh2c==1.and.ipert> natom) tim_nonlop=8
260 ! if (tim_getgh2c==2.and.ipert> natom) tim_nonlop=5
261 ! if (tim_getgh2c==3                 ) tim_nonlop=0
262 
263  idir1=alpha(idir);idir2=beta(idir)
264 
265  compute_conjugate = .false.
266  if(present(conj)) compute_conjugate = conj
267 
268 !======================================================================
269 !== Apply the 2nd-order local potential to the wavefunction
270 !======================================================================
271 
272  if (ipert/=natom+10.and.ipert/=natom+11.and.optlocal>0) then
273    msg='local part not implemented'
274    MSG_BUG(msg)
275  else
276 !  In the case of ddk operator, no local contribution (also because no self-consistency)
277 !$OMP PARALLEL DO
278    do ipw=1,npw1*my_nspinor
279      gh2c(:,ipw)=zero
280    end do
281 
282  end if
283 
284 !======================================================================
285 !== Apply the 2st-order non-local potential to the wavefunction
286 !======================================================================
287 
288  has_vnl=(ipert==natom+10.or.ipert==natom+11.or.pert_phon_elfd)
289 
290 !Use of gvnl2 depends on usevnl
291  if (usevnl==1) then
292    gvnl2_ => gvnl2
293  else
294    ABI_ALLOCATE(gvnl2_,(2,npw1*my_nspinor))
295  end if
296 
297  if (has_vnl.and.(optnl>0.or.sij_opt/=0)) then
298 
299    idirc=3*(idir1-1)+idir2 !xx=1, xy=2, xz=3, yx=4, yy=5, yz=6, zx=7, zy=8, zz=9, (xyz,xyz)=(idir1,idir2)
300 
301 ! d^2[H_nl]/dk1dk2
302 !  -------------------------------------------
303    if (ipert==natom+10) then
304      if (gs_hamkq%usepaw==1) then
305        if (usecprj==1) then
306          cwaveprj_ptr => cwaveprj
307        else
308          ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor))
309          call pawcprj_alloc(cwaveprj_tmp,0,gs_hamkq%dimcprj)
310          cwaveprj_ptr => cwaveprj_tmp
311        end if
312        cpopt=-1+5*usecprj
313        choice=8; signs=2; paw_opt=1; if (sij_opt/=0) paw_opt=sij_opt+3
314        call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/lambda/),mpi_enreg,1,nnlout,&
315 &       paw_opt,signs,gs2c,tim_nonlop,cwavef,gvnl2_)
316        if (usecprj==0) then
317          call pawcprj_free(cwaveprj_tmp)
318          ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
319        end if
320        nullify(cwaveprj_ptr)
321      else
322        choice=8; signs=2; cpopt=-1 ; paw_opt=0
323        call nonlop(choice,cpopt,cwaveprj,enlout_dum,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
324 &       paw_opt,signs,svectout_dum,tim_nonlop,cwavef,gvnl2_)
325      end if
326 
327 ! d^2[H_nl]/dk1dE2 : Non-zero only in PAW
328 !  -------------------------------------------
329    else if (ipert==natom+11.and.gs_hamkq%usepaw==1) then
330 
331      ABI_ALLOCATE(nonlop_out,(2,npw1*my_nspinor))
332 
333      if (usecprj==1) then
334        cwaveprj_ptr => cwaveprj
335      else
336        ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor))
337        call pawcprj_alloc(cwaveprj_tmp,2,gs_hamkq%dimcprj)
338        cwaveprj_ptr => cwaveprj_tmp
339      end if
340 
341      if (opt_gvnl2==1.and.optnl>=1) then
342 
343 !      Compute application of dS/dk1 to i*d[cwavef]/dk2
344 !      sum_{i,j} s_ij d(|p_i><p_j|)/dk(idir1) | i*psi^(k(idir2)) >
345        cpopt=-1 ; choice=5 ; paw_opt=3 ; signs=2
346        call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
347 &       paw_opt,signs,nonlop_out,tim_nonlop,gvnl2_,vectout_dum)
348 
349 !$OMP PARALLEL DO
350        do ipw=1,npw1*my_nspinor
351          gvnl2_(:,ipw)=nonlop_out(:,ipw)
352        end do
353 
354 !      Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
355 !      sum_{i,j} chi_ij(idir2) d(|p_i><p_j|)/dk(idir1) | psi^(0) >
356        cpopt=4*usecprj ; choice=5 ; paw_opt=1 ; signs=2
357        call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
358 &       paw_opt,signs,svectout_dum,tim_nonlop,cwavef,nonlop_out,enl=enl_ptr)
359 
360 !$OMP PARALLEL DO
361        do ipw=1,npw1*my_nspinor
362          gvnl2_(:,ipw)=gvnl2_(:,ipw)+nonlop_out(:,ipw)
363        end do
364 
365      else
366 
367 !$OMP PARALLEL DO
368        do ipw=1,npw1*my_nspinor
369          gvnl2_(:,ipw)=zero
370        end do
371 
372      end if ! opt_gvnl2==1
373 
374 !    Compute derivatives due to projectors |d^2[p_i]/dk1dk2>,|d[p_i]/dk1>,|d[p_i]/dk2>
375 !    i * sum_{i,j} (d(|p_i><dp_j/dk(idir2)|)/dk(idir1) | psi^(0) >
376      cpopt=-1+5*usecprj ; choice=81 ; paw_opt=3 ; signs=2
377      call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
378 &     paw_opt,signs,nonlop_out,tim_nonlop,cwavef,vectout_dum)
379 
380      if(compute_conjugate) then
381 !$OMP PARALLEL DO
382        do ipw=1,npw1*my_nspinor ! Note the multiplication by -i
383          gvnl2_(1,ipw)=gvnl2_(1,ipw)+nonlop_out(2,ipw)
384          gvnl2_(2,ipw)=gvnl2_(2,ipw)-nonlop_out(1,ipw)
385        end do
386      else
387 !$OMP PARALLEL DO
388        do ipw=1,npw1*my_nspinor ! Note the multiplication by i
389          gvnl2_(1,ipw)=gvnl2_(1,ipw)-nonlop_out(2,ipw)
390          gvnl2_(2,ipw)=gvnl2_(2,ipw)+nonlop_out(1,ipw)
391        end do
392      end if
393 
394      ABI_DEALLOCATE(nonlop_out)
395      if (sij_opt==1) gs2c=zero
396      if (usecprj==0) then
397        call pawcprj_free(cwaveprj_tmp)
398        ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
399      end if
400      nullify(cwaveprj_ptr)
401 
402 ! d^2[H_nl]/dtau1dE2 : Non-zero only in PAW
403 !  -------------------------------------------
404    else if (pert_phon_elfd.and.gs_hamkq%usepaw==1) then
405 
406      iatm = ipert-(natom+11)
407      if (iatm<1.or.iatm>natom) then
408        MSG_BUG(" iatm must be between 1 and natom")
409      end if
410 
411      ABI_ALLOCATE(nonlop_out,(2,npw1*my_nspinor))
412 
413      if (usecprj==1) then
414        cwaveprj_ptr => cwaveprj
415      else
416        ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor))
417        call pawcprj_alloc(cwaveprj_tmp,2,gs_hamkq%dimcprj)
418        cwaveprj_ptr => cwaveprj_tmp
419      end if
420 
421      if (opt_gvnl2==1) then
422 
423 !      Compute application of dS/dtau1 to i*d[cwavef]/dk2
424 !      sum_{i,j} s_ij d(|p_i><p_j|)/dtau(idir1) | i*psi^(k(idir2)) >
425        cpopt=-1 ; choice=2 ; paw_opt=3 ; signs=2
426        call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
427 &       paw_opt,signs,nonlop_out,tim_nonlop,gvnl2_,vectout_dum,iatom_only=iatm)
428 
429 !$OMP PARALLEL DO
430        do ipw=1,npw1*my_nspinor
431          gvnl2_(:,ipw)=nonlop_out(:,ipw)
432        end do
433 
434 !      Compute part of H^(2) due to derivative of projectors (idir1) and derivative of Dij (idir2)
435 !      sum_{i,j} chi_ij(idir2) d(|p_i><p_j|)/dtau(idir1) | psi^(0) >
436        cpopt=4*usecprj ; choice=2 ; paw_opt=1 ; signs=2
437        call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idir1,(/zero/),mpi_enreg,1,nnlout,&
438 &       paw_opt,signs,svectout_dum,tim_nonlop,cwavef,nonlop_out,enl=enl_ptr,iatom_only=iatm)
439 
440 !$OMP PARALLEL DO
441        do ipw=1,npw1*my_nspinor
442          gvnl2_(:,ipw)=gvnl2_(:,ipw)+nonlop_out(:,ipw)
443        end do
444 
445      else
446 
447 !$OMP PARALLEL DO
448        do ipw=1,npw1*my_nspinor
449          gvnl2_(:,ipw)=zero
450        end do
451 
452      end if ! opt_gvnl2==1
453 
454 !    Compute derivatives due to projectors |d^2[p_i]/dtau1dk2>,|d[p_i]/dtau1>,|d[p_i]/dk2>
455 !    i * sum_{i,j} (d(|p_i><dp_j/dk(idir2)|)/dtau(idir1) | psi^(0) >
456      cpopt=-1+5*usecprj ; choice=54 ; paw_opt=3 ; signs=2
457      call nonlop(choice,cpopt,cwaveprj_ptr,enlout_dum,gs_hamkq,idirc,(/zero/),mpi_enreg,1,nnlout,&
458 &     paw_opt,signs,nonlop_out,tim_nonlop,cwavef,vectout_dum,iatom_only=iatm)
459 
460      if(compute_conjugate) then
461 !$OMP PARALLEL DO
462        do ipw=1,npw1*my_nspinor ! Note the multiplication by -i
463          gvnl2_(1,ipw)=gvnl2_(1,ipw)+nonlop_out(2,ipw)
464          gvnl2_(2,ipw)=gvnl2_(2,ipw)-nonlop_out(1,ipw)
465        end do
466      else
467 !$OMP PARALLEL DO
468        do ipw=1,npw1*my_nspinor ! Note the multiplication by i
469          gvnl2_(1,ipw)=gvnl2_(1,ipw)-nonlop_out(2,ipw)
470          gvnl2_(2,ipw)=gvnl2_(2,ipw)+nonlop_out(1,ipw)
471        end do
472      end if
473 
474      ABI_DEALLOCATE(nonlop_out)
475      if (sij_opt==1) gs2c=zero
476      if (usecprj==0) then
477        call pawcprj_free(cwaveprj_tmp)
478        ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
479      end if
480      nullify(cwaveprj_ptr)
481 
482    end if
483 
484 !No non-local part
485 !-------------------------------------------
486  else
487 
488    if (optnl>=1) then
489  !$OMP PARALLEL DO
490      do ipw=1,npw1*my_nspinor
491        gvnl2_(:,ipw)=zero
492      end do
493    end if
494    if (sij_opt/=0) then
495  !$OMP PARALLEL DO
496      do ipw=1,npw1*my_nspinor
497        gs2c(:,ipw)=zero
498      end do
499    end if
500 
501  end if
502 
503  if (associated(enl_ptr)) then
504    nullify(enl_ptr)
505  end if
506  if (allocated(enl_temp)) then
507    ABI_DEALLOCATE(enl_temp)
508  end if
509 
510 !======================================================================
511 !== Apply the 2nd-order kinetic operator to the wavefunction
512 !======================================================================
513 
514  if (present(optkin)) then
515    has_kin=(optkin/=0.and.ipert==natom+10)
516  else
517    has_kin=(ipert==natom+10)
518  end if
519 
520 !k-point perturbation
521 !-------------------------------------------
522  if (associated(gs_hamkq%kinpw_kp)) then
523    kinpw1 => gs_hamkq%kinpw_kp
524  else if (optnl>=1.or.has_kin) then
525    msg='need kinpw1 allocated!'
526    MSG_BUG(msg)
527  end if
528  if (associated(rf_hamkq%ddkinpw_k)) then
529    ddkinpw => rf_hamkq%ddkinpw_k
530  else if (has_kin) then
531    msg='need ddkinpw allocated!'
532    MSG_BUG(msg)
533  end if
534 
535  if (has_kin) then
536    do ispinor=1,my_nspinor
537  !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwavef,ispinor,gvnl2_,ddkinpw,kinpw1,npw,my_nspinor)
538      do ipw=1,npw
539        ipws=ipw+npw*(ispinor-1)
540        if(kinpw1(ipw)<huge(zero)*1.d-11)then
541          gvnl2_(1,ipws)=gvnl2_(1,ipws)+ddkinpw(ipw)*cwavef(1,ipws)
542          gvnl2_(2,ipws)=gvnl2_(2,ipws)+ddkinpw(ipw)*cwavef(2,ipws)
543        else
544          gvnl2_(1,ipws)=zero
545          gvnl2_(2,ipws)=zero
546        end if
547      end do
548    end do
549  end if
550 
551 !======================================================================
552 !== Sum contributions to get the application of H^(2) to the wf
553 !======================================================================
554 !Also filter the wavefunctions for large modified kinetic energy
555 
556 !Add non-local+kinetic to local part
557  if (optnl>=1.or.has_kin) then
558    do ispinor=1,my_nspinor
559      ipws=(ispinor-1)*npw1
560  !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh2c,gvnl2_,kinpw1,ipws,npw1)
561      do ipw=1+ipws,npw1+ipws
562        if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then
563          gh2c(1,ipw)=gh2c(1,ipw)+gvnl2_(1,ipw)
564          gh2c(2,ipw)=gh2c(2,ipw)+gvnl2_(2,ipw)
565        else
566          gh2c(1,ipw)=zero
567          gh2c(2,ipw)=zero
568        end if
569      end do
570    end do
571  end if
572 
573  if (usevnl==1) then
574    nullify(gvnl2_)
575  else
576    ABI_DEALLOCATE(gvnl2_)
577  end if
578 
579 !call timab(196+tim_getgh2c,2,tsec)
580 
581  DBG_EXIT("COLL")
582 
583 end subroutine getgh2c

ABINIT/m_getgh2c [ Modules ]

[ Top ] [ Modules ]

NAME

  m_getgh2c

FUNCTION

COPYRIGHT

  Copyright (C) 2015-2018 ABINIT group (MT,JLJ)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_getgh2c
28 
29  use defs_basis
30  use defs_abitypes
31  use m_abicore
32  use m_errors
33 
34  use m_pawcprj,     only : pawcprj_type,pawcprj_alloc,pawcprj_free
35  use m_hamiltonian, only : gs_hamiltonian_type,rf_hamiltonian_type
36  use m_nonlop,      only : nonlop
37 
38  implicit none
39 
40  private