TABLE OF CONTENTS


ABINIT/getdc1 [ Functions ]

[ Top ] [ Functions ]

NAME

 getdc1

FUNCTION

 Compute |delta_C^(1)> from one wave function C - PAW ONLY
 Compute <G|delta_C^(1)> and eventually <P_i| delta_C^(1)> (P_i= non-local projector)
 delta_C^(1) is the variation of wavefunction only due to variation of overlap operator S.
 delta_C^(1)=-1/2.Sum_j [ <C_j|S^(1)|C>.C_j
         see PRB 78, 035105 (2008) [[cite:Audouze2008]], Eq. (42)

INPUTS

  cgq(2,mcgq)=wavefunction coefficients for ALL bands at k+Q
  cprjq(natom,mcprjq)= wave functions at k+q projected with non-local projectors: cprjq=<P_i|Cnk+q>
  ibgq=shift to be applied on the location of data in the array cprjq
  icgq=shift to be applied on the location of data in the array cgq
  istwfk=option parameter that describes the storage of wfs
  mcgq=second dimension of the cgq array
  mcprjq=second dimension of the cprjq array
  mpi_enreg=information about MPI parallelization
  natom= number of atoms in cell
  nband=number of bands
  npw1=number of planewaves in basis sphere at k+Q
  nspinor=number of spinorial components of the wavefunctions
  opt_cprj=flag governing the computation of <P_i|delta_C^(1)> (P_i= non-local projector)
  s1cwave0(2,npw1*nspinor)=<G|S^(1)|C> where S^(1) is the first-order overlap operator

OUTPUT

  dcwavef(2,npw1*nspinor)=change of wavefunction due to change of overlap PROJECTED ON PLANE-WAVES:
         dcwavef is delta_C(1)=-1/2.Sum_{j}[<C0_k+q_j|S(1)|C0_k_i>.|C0_k+q_j>]
  === if optcprj=1 ===
  dcwaveprj(natom,nspinor*optcprj)=change of wavefunction due to change of overlap PROJECTED ON NL-PROJECTORS:

PARENTS

      dfpt_cgwf,dfpt_nstpaw

CHILDREN

      pawcprj_axpby,pawcprj_lincom,projbd

SOURCE

1254 subroutine getdc1(cgq,cprjq,dcwavef,dcwaveprj,ibgq,icgq,istwfk,mcgq,mcprjq,&
1255 &                 mpi_enreg,natom,nband,npw1,nspinor,optcprj,s1cwave0)
1256 
1257 
1258 !This section has been created automatically by the script Abilint (TD).
1259 !Do not modify the following lines by hand.
1260 #undef ABI_FUNC
1261 #define ABI_FUNC 'getdc1'
1262 !End of the abilint section
1263 
1264  implicit none
1265 
1266 !Arguments ------------------------------------
1267 !scalars
1268  integer,intent(in) :: ibgq,icgq,istwfk,mcgq,mcprjq,natom,nband,npw1,nspinor,optcprj
1269  type(MPI_type),intent(in) :: mpi_enreg
1270 !arrays
1271  real(dp),intent(in) :: cgq(2,mcgq),s1cwave0(2,npw1*nspinor)
1272  real(dp),intent(out) :: dcwavef(2,npw1*nspinor)
1273  type(pawcprj_type),intent(in) :: cprjq(natom,mcprjq)
1274  type(pawcprj_type),intent(inout) :: dcwaveprj(natom,nspinor*optcprj)
1275 
1276 !Local variables-------------------------------
1277 !scalars
1278  integer, parameter :: tim_projbd=0
1279  integer :: ipw
1280  real(dp),parameter :: scal=-half
1281 !arrays
1282  real(dp), allocatable :: dummy(:,:),scprod(:,:)
1283  type(pawcprj_type),allocatable :: tmpcprj(:,:)
1284 
1285 ! *********************************************************************
1286 
1287  DBG_ENTER("COLL")
1288 
1289 !$OMP PARALLEL DO
1290  do ipw=1,npw1*nspinor
1291    dcwavef(1:2,ipw)=s1cwave0(1:2,ipw)
1292  end do
1293 
1294  ABI_ALLOCATE(dummy,(0,0))
1295  ABI_ALLOCATE(scprod,(2,nband))
1296 
1297 !=== 1- COMPUTE: <G|S^(1)|C_k> - Sum_j [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>]
1298 !!               using the projb routine
1299 !Note the subtlety: projbd is called with useoverlap=0 and s1cwave0
1300 !in order to get Sum[<cgq|s1|c>|cgq>]=Sum[<cgq|gs1>|cgq>]
1301  call projbd(cgq,dcwavef,-1,icgq,0,istwfk,mcgq,0,nband,npw1,nspinor,&
1302 & dummy,scprod,0,tim_projbd,0,mpi_enreg%me_g0,mpi_enreg%comm_fft)
1303 
1304 !=== 2- COMPUTE: <G|delta_C^(1)> = -1/2.Sum_j [<C_k+q,j|S^(1)|C_k>.<G|C_k+q,j>] by substraction
1305 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(dcwavef,s1cwave0,npw1,nspinor)
1306  do ipw=1,npw1*nspinor
1307    dcwavef(1:2,ipw)=scal*(s1cwave0(1:2,ipw)-dcwavef(1:2,ipw))
1308  end do
1309 
1310 !=== 3- COMPUTE: <P_i|delta_C^(1)> = -1/2.Sum_j [<C_k+q,j|S^(1)|C_k>.<P_i|C_k+q,j>]
1311  if (optcprj==1.and.mcprjq>0) then
1312    ABI_DATATYPE_ALLOCATE(tmpcprj,(natom,nspinor))
1313    call pawcprj_lincom(scprod,cprjq(:,ibgq+1:ibgq+nspinor*nband),dcwaveprj,nband)
1314    call pawcprj_axpby(zero,scal,tmpcprj,dcwaveprj)
1315    ABI_DATATYPE_DEALLOCATE(tmpcprj)
1316  end if
1317 
1318  ABI_DEALLOCATE(dummy)
1319  ABI_DEALLOCATE(scprod)
1320 
1321  DBG_EXIT("COLL")
1322 
1323 end subroutine getdc1

ABINIT/getgh1c [ Functions ]

[ Top ] [ Functions ]

NAME

 getgh1c

FUNCTION

 Compute <G|H^(1)|C> (or <G|H^(1)-lambda.S^(1)|C>) for input vector |C> expressed in reciprocal space.
 (H^(1) is the 1st-order pertubed Hamiltonian, S^(1) is the 1st-order perturbed overlap operator).
 Result is put in array gh1c.
 If required, part of <G|K(1)+Vnonlocal^(1)|C> not depending on VHxc^(1) is also returned in gvnl1c.
 If required, <G|S^(1)|C> is returned in gs1c (S=overlap - PAW only)

INPUTS

  berryopt=option for Berry phase
  cwave(2,npw*nspinor)=input wavefunction, in reciprocal space
  cwaveprj(natom,nspinor*usecprj)=<p_lmn|C> coefficients for wavefunction |C> (and 1st derivatives)
     if not allocated or size=0, they are locally computed (and not sotred)
  dkinpw(npw)=derivative of the (modified) kinetic energy for each plane wave at k (Hartree)
  grad_berry(2,npw1*nspinor*(berryopt/4))= the gradient of the Berry phase term
  gs_hamkq <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k+q
  idir=direction of the perturbation
  ipert=type of the perturbation
  lambda=real use to apply H^(1)-lambda.S^(1)
  mpi_enreg=information about MPI parallelization
  npw=number of planewaves in basis sphere at given k.
  npw1=number of planewaves in basis sphere at k+q
  optlocal=0: local part of H^(1) is not computed in gh1c=<G|H^(1)|C>
           1: local part of H^(1) is computed in gh1c=<G|H^(1)|C>
  optnl=0: non-local part of H^(1) is not computed in gh1c=<G|H^(1)|C>
        1: non-local part of H^(1) depending on VHxc^(1) is not computed in gh1c=<G|H^(1)|C>
        2: non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C>
  opt_gvnl1=option controlling the use of gvnl1 array:
            0: used as an output
            1: used as an input:   (only for ipert=natom+2)
                 NCPP: contains the ddk 1-st order WF
                 PAW: contains frozen part of 1st-order hamiltonian
            2: used as input/ouput:    - used only for PAW and ipert=natom+2
                 At input: contains the ddk 1-st order WF (times i)
                 At output: contains frozen part of 1st-order hamiltonian
  rf_hamkq <type(rf_hamiltonian_type)>=all data for the 1st-order Hamiltonian at k,k+q
  sij_opt= -PAW ONLY-  if  0, only matrix elements <G|H^(1)|C> have to be computed
     (S=overlap)       if  1, matrix elements <G|S^(1)|C> have to be computed in gs1c in addition to gh1c
                       if -1, matrix elements <G|H^(1)-lambda.S^(1)|C> have to be computed in gh1c (gs1c not used)
  tim_getgh1c=timing code of the calling subroutine (can be set to 0 if not attributed)
  usevnl=1 if gvnl1=(part of <G|K^(1)+Vnl^(1)-lambda.S^(1)|C> not depending on VHxc^(1)) has to be input/output

OUTPUT

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

PARENTS

      dfpt_cgwf,dfpt_nstpaw,dfpt_nstwf,dfpt_wfkfermi,m_gkk,m_phgamma,m_phpi
      m_rf2,m_sigmaph

CHILDREN

      kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian,load_kprime_hamiltonian
      mkffnl,mkkin,mkkpg

SOURCE

127 subroutine getgh1c(berryopt,cwave,cwaveprj,gh1c,grad_berry,gs1c,gs_hamkq,&
128 &          gvnl1,idir,ipert,lambda,mpi_enreg,optlocal,optnl,opt_gvnl1,&
129 &          rf_hamkq,sij_opt,tim_getgh1c,usevnl,conj)
130 
131 
132 !This section has been created automatically by the script Abilint (TD).
133 !Do not modify the following lines by hand.
134 #undef ABI_FUNC
135 #define ABI_FUNC 'getgh1c'
136 !End of the abilint section
137 
138  implicit none
139 
140 !Arguments ------------------------------------
141 !scalars
142  logical,intent(in),optional :: conj
143  integer,intent(in) :: berryopt,idir,ipert,optlocal,optnl,opt_gvnl1,sij_opt,tim_getgh1c,usevnl
144  real(dp),intent(in) :: lambda
145  type(MPI_type),intent(in) :: mpi_enreg
146  type(gs_hamiltonian_type),intent(inout),target :: gs_hamkq
147  type(rf_hamiltonian_type),intent(inout),target :: rf_hamkq
148 !arrays
149  real(dp),intent(in) :: grad_berry(:,:)
150  real(dp),intent(inout) :: cwave(2,gs_hamkq%npw_k*gs_hamkq%nspinor)
151  real(dp),intent(out) :: gh1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor)
152  real(dp),intent(out) :: gs1c(2,gs_hamkq%npw_kp*gs_hamkq%nspinor)
153  real(dp),intent(inout),target :: gvnl1(2,gs_hamkq%npw_kp*gs_hamkq%nspinor)
154  type(pawcprj_type),intent(inout),target :: cwaveprj(:,:)
155 
156 !Local variables-------------------------------
157 !scalars
158  integer,parameter :: level=16
159  integer :: choice,cplex1,cpopt,ipw,ipws,ispinor,istr,i1,i2,i3
160  integer :: my_nspinor,natom,ncpgr,nnlout=1,npw,npw1,paw_opt,signs
161  integer :: tim_fourwf,tim_nonlop,usecprj
162  logical :: compute_conjugate,has_kin,usevnl2
163  real(dp) :: weight
164  character(len=500) :: msg
165 !arrays
166  real(dp) :: enlout(1),tsec(2),svectout_dum(1,1),vectout_dum(1,1)
167  real(dp),allocatable :: cwave_sp(:,:),cwavef1(:,:),cwavef2(:,:)
168  real(dp),allocatable :: gh1c_sp(:,:),gh1c1(:,:),gh1c2(:,:),gh1c3(:,:),gh1c4(:,:),gvnl2(:,:)
169  real(dp),allocatable :: nonlop_out(:,:),vlocal1_tmp(:,:,:),work(:,:,:,:)
170  real(dp),ABI_CONTIGUOUS pointer :: gvnl1_(:,:)
171  real(dp),pointer :: dkinpw(:),kinpw1(:)
172  type(pawcprj_type),allocatable,target :: cwaveprj_tmp(:,:)
173  type(pawcprj_type),pointer :: cwaveprj_ptr(:,:)
174 
175 ! *********************************************************************
176 
177  DBG_ENTER("COLL")
178 
179 !Keep track of total time spent in getgh1c
180  call timab(196+tim_getgh1c,1,tsec)
181 
182 !======================================================================
183 !== Initialisations and compatibility tests
184 !======================================================================
185 
186  npw  =gs_hamkq%npw_k
187  npw1 =gs_hamkq%npw_kp
188  natom=gs_hamkq%natom
189 
190 !Compatibility tests
191  if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then
192    if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))).or. &
193 &      (optnl==2.and.(.not.associated(rf_hamkq%e1kbsc))))then
194      msg='ekb derivatives must be allocated for ipert<=natom or natom+3/4 !'
195      MSG_BUG(msg)
196    end if
197  end if
198  if(gs_hamkq%usepaw==1.and.(ipert==natom+2)) then
199    if ((optnl>=1.and.(.not.associated(rf_hamkq%e1kbfr))).or. &
200 &      (optnl==2.and.(.not.associated(rf_hamkq%e1kbsc))))then
201      msg='ekb derivatives must be allocated for ipert=natom+2 !'
202      MSG_BUG(msg)
203    end if
204    if (usevnl==0) then
205      msg='gvnl1 must be allocated for ipert=natom+2 !'
206      MSG_BUG(msg)
207    end if
208  end if
209  if(ipert==natom+2.and.opt_gvnl1==0) then
210    msg='opt_gvnl1=0 not compatible with ipert=natom+2 !'
211    MSG_BUG(msg)
212  end if
213  if (mpi_enreg%paral_spinor==1) then
214    msg='Not compatible with parallelization over spinorial components !'
215    MSG_BUG(msg)
216  end if
217 
218 !Check sizes
219  my_nspinor=max(1,gs_hamkq%nspinor/mpi_enreg%nproc_spinor)
220  if (size(cwave)<2*npw*my_nspinor) then
221    msg='wrong size for cwave!'
222    MSG_BUG(msg)
223  end if
224  if (size(gh1c)<2*npw1*my_nspinor) then
225    msg='wrong size for gh1c!'
226    MSG_BUG(msg)
227  end if
228  if (usevnl/=0) then
229    if (size(gvnl1)<2*npw1*my_nspinor) then
230      msg='wrong size for gvnl1!'
231      MSG_BUG(msg)
232    end if
233  end if
234  if (sij_opt==1) then
235    if (size(gs1c)<2*npw1*my_nspinor) then
236      msg='wrong size for gs1c!'
237      MSG_BUG(msg)
238    end if
239  end if
240  if (berryopt>=4) then
241    if (size(grad_berry)<2*npw1*my_nspinor) then
242      msg='wrong size for grad_berry!'
243      MSG_BUG(msg)
244    end if
245  end if
246 
247 !PAW: specific treatment for usecprj input arg
248 !     force it to zero if cwaveprj is not allocated
249  usecprj=gs_hamkq%usecprj ; ncpgr=0
250  if(gs_hamkq%usepaw==1) then
251    if (size(cwaveprj)==0) usecprj=0
252    if (usecprj/=0) then
253      ncpgr=cwaveprj(1,1)%ncpgr
254      if (size(cwaveprj)<gs_hamkq%natom*my_nspinor) then
255        msg='wrong size for cwaveprj!'
256        MSG_BUG(msg)
257      end if
258      if(gs_hamkq%usepaw==1.and.(ipert>=0.and.(ipert<=natom.or.ipert==natom+3.or.ipert==natom+4))) then
259        if (ncpgr/=1)then
260          msg='Projected WFs (cprj) derivatives are not correctly stored !'
261          MSG_BUG(msg)
262        end if
263      end if
264    end if
265  else
266    if(usecprj==1)then
267      msg='usecprj==1 not allowed for NC psps !'
268      MSG_BUG(msg)
269    end if
270  end if
271 
272  tim_nonlop=8
273  if (tim_getgh1c==1.and.ipert<=natom) tim_nonlop=7
274  if (tim_getgh1c==2.and.ipert<=natom) tim_nonlop=5
275  if (tim_getgh1c==1.and.ipert> natom) tim_nonlop=8
276  if (tim_getgh1c==2.and.ipert> natom) tim_nonlop=5
277  if (tim_getgh1c==3                 ) tim_nonlop=0
278 
279  compute_conjugate = .false.
280  if(present(conj)) compute_conjugate = conj
281 
282 !======================================================================
283 !== Apply the 1st-order local potential to the wavefunction
284 !======================================================================
285 
286 !Phonon perturbation
287 !or Electric field perturbation
288 !or Strain perturbation
289 !-------------------------------------------
290  if (ipert<=natom+5.and.ipert/=natom+1.and.optlocal>0) then !SPr deb
291 
292    ABI_ALLOCATE(work,(2,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
293 
294    if (gs_hamkq%nvloc==1) then
295 
296      weight=one ; tim_fourwf=4
297      call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave,gh1c,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
298 &     gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
299 &     npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
300 &     use_gpu_cuda=gs_hamkq%use_gpu_cuda)
301      if(gs_hamkq%nspinor==2)then
302        ABI_ALLOCATE(cwave_sp,(2,npw))
303        ABI_ALLOCATE(gh1c_sp,(2,npw1))
304 !$OMP PARALLEL DO
305        do ipw=1,npw
306          cwave_sp(1,ipw)=cwave(1,ipw+npw)
307          cwave_sp(2,ipw)=cwave(2,ipw+npw)
308        end do
309        call fourwf(rf_hamkq%cplex,rf_hamkq%vlocal1,cwave_sp,gh1c_sp,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
310 &       gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
311 &       npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
312 &       use_gpu_cuda=gs_hamkq%use_gpu_cuda)
313 !$OMP PARALLEL DO
314        do ipw=1,npw1
315          gh1c(1,ipw+npw1)=gh1c_sp(1,ipw)
316          gh1c(2,ipw+npw1)=gh1c_sp(2,ipw)
317        end do
318        ABI_DEALLOCATE(cwave_sp)
319        ABI_DEALLOCATE(gh1c_sp)
320      end if
321    else ! Non-Collinear magnetism for nvloc=4
322      if (gs_hamkq%nspinor==2) then
323        weight=one ; tim_fourwf=4
324        ABI_ALLOCATE(gh1c1,(2,npw1))
325        ABI_ALLOCATE(gh1c2,(2,npw1))
326        ABI_ALLOCATE(gh1c3,(2,npw1))
327        ABI_ALLOCATE(gh1c4,(2,npw1))
328        gh1c1(:,:)=zero; gh1c2(:,:)=zero; gh1c3(:,:)=zero ;  gh1c4(:,:)=zero
329        ABI_ALLOCATE(vlocal1_tmp,(rf_hamkq%cplex*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6)) !SPr: notation/dimension corrected vlocal_tmp -> vlocal1_tmp
330        ABI_ALLOCATE(cwavef1,(2,npw))
331        ABI_ALLOCATE(cwavef2,(2,npw))
332        do ipw=1,npw
333          cwavef1(1:2,ipw)=cwave(1:2,ipw)
334          cwavef2(1:2,ipw)=cwave(1:2,ipw+npw)
335        end do
336 !      gh1c1=v11*phi1
337        vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,1)
338        call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef1,gh1c1,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
339 &       gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
340 &       npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
341 &       use_gpu_cuda=gs_hamkq%use_gpu_cuda)
342 !      gh1c2=v22*phi2
343        vlocal1_tmp(:,:,:)=rf_hamkq%vlocal1(:,:,:,2)
344        call fourwf(rf_hamkq%cplex,vlocal1_tmp,cwavef2,gh1c2,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
345 &       gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
346 &       npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
347 &       use_gpu_cuda=gs_hamkq%use_gpu_cuda)
348        ABI_DEALLOCATE(vlocal1_tmp)
349        cplex1=2
350        ABI_ALLOCATE(vlocal1_tmp,(cplex1*gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6))
351 !      gh1c3=(re(v12)-im(v12))*phi1 => v^21*phi1
352        if(rf_hamkq%cplex==1) then
353          do i3=1,gs_hamkq%n6
354            do i2=1,gs_hamkq%n5
355              do i1=1,gs_hamkq%n4
356                vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(i1,i2,i3,3)
357                vlocal1_tmp(2*i1  ,i2,i3)=-rf_hamkq%vlocal1(i1,i2,i3,4)
358              end do
359            end do
360          end do
361        else
362        !SPr: modified definition of local potential components for cplex=2 (see dotprod_vn)
363        !also, v21==v12* not always holds (e.g. magnetic field perturbation)
364          do i3=1,gs_hamkq%n6
365            do i2=1,gs_hamkq%n5
366              do i1=1,gs_hamkq%n4
367                vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(2*i1  ,i2,i3,4)
368                vlocal1_tmp(2*i1  ,i2,i3)=-rf_hamkq%vlocal1(2*i1-1,i2,i3,4)
369              end do
370            end do
371          end do
372        end if
373        call fourwf(cplex1,vlocal1_tmp,cwavef1,gh1c3,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
374 &       gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
375 &       npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
376 &       use_gpu_cuda=gs_hamkq%use_gpu_cuda)
377 !      gh1c4=(re(v12)+im(v12))*phi2 => v^12*phi2
378        if(rf_hamkq%cplex==1) then
379          do i3=1,gs_hamkq%n6
380            do i2=1,gs_hamkq%n5
381              do i1=1,gs_hamkq%n4
382                vlocal1_tmp(2*i1,i2,i3)=-vlocal1_tmp(2*i1,i2,i3)
383              end do
384            end do
385          end do
386        else
387          !for cplex=2 and time-reversal breaking perturbations,v21/=v12*
388          do i3=1,gs_hamkq%n6
389            do i2=1,gs_hamkq%n5
390              do i1=1,gs_hamkq%n4
391                vlocal1_tmp(2*i1-1,i2,i3)= rf_hamkq%vlocal1(2*i1-1,i2,i3,3)
392                vlocal1_tmp(2*i1  ,i2,i3)= rf_hamkq%vlocal1(2*i1  ,i2,i3,3)
393              end do
394            end do
395          end do
396        end if
397        call fourwf(cplex1,vlocal1_tmp,cwavef2,gh1c4,work,gs_hamkq%gbound_k,gs_hamkq%gbound_kp,&
398 &       gs_hamkq%istwf_k,gs_hamkq%kg_k,gs_hamkq%kg_kp,gs_hamkq%mgfft,mpi_enreg,1,gs_hamkq%ngfft,&
399 &       npw,npw1,gs_hamkq%n4,gs_hamkq%n5,gs_hamkq%n6,2,mpi_enreg%paral_kgb,tim_fourwf,weight,weight,&
400 &       use_gpu_cuda=gs_hamkq%use_gpu_cuda)
401        ABI_DEALLOCATE(vlocal1_tmp)
402 !      Build gh1c from pieces
403 !      gh1c_1 = (v11, v12) (psi1) matrix vector product
404 !      gh1c_2 = (v12*,v22) (psi2)
405        do ipw=1,npw1
406          gh1c(1:2,ipw)     =gh1c1(1:2,ipw)+gh1c4(1:2,ipw)
407          gh1c(1:2,ipw+npw1)=gh1c3(1:2,ipw)+gh1c2(1:2,ipw)
408        end do
409        ABI_DEALLOCATE(gh1c1)
410        ABI_DEALLOCATE(gh1c2)
411        ABI_DEALLOCATE(gh1c3)
412        ABI_DEALLOCATE(gh1c4)
413        ABI_DEALLOCATE(cwavef1)
414        ABI_DEALLOCATE(cwavef2)
415      else
416        msg='nspinor/=1 for Non-collinear calculations!'
417        MSG_BUG(msg)
418      end if
419    end if ! nvloc
420 
421    ABI_DEALLOCATE(work)
422 
423 !  k-point perturbation (or no local part, i.e. optlocal=0)
424 !  -------------------------------------------
425  else if (ipert==natom+1.or.optlocal==0) then
426 
427 !  In the case of ddk operator, no local contribution (also because no self-consistency)
428 !$OMP PARALLEL DO
429    do ipw=1,npw1*my_nspinor
430      gh1c(:,ipw)=zero
431    end do
432 
433  end if
434 
435 !======================================================================
436 !== Apply the 1st-order non-local potential to the wavefunction
437 !======================================================================
438 
439 !Use of gvnl1 depends on usevnl
440  if (usevnl==1) then
441    gvnl1_ => gvnl1
442  else
443    ABI_ALLOCATE(gvnl1_,(2,npw1*my_nspinor))
444  end if
445 
446 !Phonon perturbation
447 !-------------------------------------------
448  if (ipert<=natom.and.(optnl>0.or.sij_opt/=0)) then
449 
450 !  PAW:
451    if (gs_hamkq%usepaw==1) then
452 
453      if (usecprj==1) then
454        cwaveprj_ptr => cwaveprj
455      else
456        ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor))
457        call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
458        cwaveprj_ptr => cwaveprj_tmp
459      end if
460 
461 !    1- Compute derivatives due to projectors |p_i>^(1)
462 !    Only displaced atom contributes
463      cpopt=-1+5*usecprj ; choice=2 ; signs=2
464      paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3
465      call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
466 &     paw_opt,signs,gs1c,tim_nonlop,cwave,gvnl1_,iatom_only=ipert)
467 
468 !    2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1))
469 !    All atoms contribute
470      if (optnl>=1) then
471        ABI_ALLOCATE(nonlop_out,(2,npw1*my_nspinor))
472        cpopt=1+3*usecprj ; choice=1 ; signs=2 ; paw_opt=1
473        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
474 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr)
475 !$OMP PARALLEL DO
476        do ipw=1,npw1*my_nspinor
477          gvnl1_(:,ipw)=gvnl1_(:,ipw)+nonlop_out(:,ipw)
478        end do
479        ABI_DEALLOCATE(nonlop_out)
480      end if
481 
482 !    3- Compute derivatives due to self-consistent part of D_ij^(1) (depending on VHxc^(1))
483 !    All atoms contribute
484      if (optnl==2) then
485        ABI_ALLOCATE(gvnl2,(2,npw1*my_nspinor))
486        cpopt=4 ; choice=1 ; signs=2 ; paw_opt=1
487        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
488 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc)
489      end if
490 
491      if (usecprj==0) then
492        call pawcprj_free(cwaveprj_tmp)
493        ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
494      end if
495      nullify(cwaveprj_ptr)
496 
497 !  Norm-conserving psps:
498    else
499 !    Compute only derivatives due to projectors |p_i>^(1)
500      cpopt=-1 ; choice=2 ; signs=2 ; paw_opt=0
501      call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
502 &     paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1_,iatom_only=ipert)
503      if (sij_opt==1) then
504 !$OMP PARALLEL DO
505        do ipw=1,npw1*my_nspinor
506          gs1c(:,ipw)=zero
507        end do
508      end if
509    end if
510 
511 !  k-point perturbation
512 !  -------------------------------------------
513  else if (ipert==natom+1.and.(optnl>0.or.sij_opt/=0)) then
514 
515    tim_nonlop=8 ; signs=2 ; choice=5
516    if (gs_hamkq%usepaw==1) then
517      if (usecprj==1) then
518        cwaveprj_ptr => cwaveprj
519      else
520        ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor))
521        call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
522        cwaveprj_ptr => cwaveprj_tmp
523      end if
524      cpopt=-1+5*usecprj; paw_opt=1; if (sij_opt/=0) paw_opt=sij_opt+3
525 !    JLJ: BUG (wrong result) of H^(1) if stored cprj are used in PAW DDKs with nspinor==2 (==1 works fine).
526 !    To be debugged, if someone has time...
527      if(gs_hamkq%nspinor==2) cpopt=-1
528      call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
529 &     paw_opt,signs,gs1c,tim_nonlop,cwave,gvnl1_)
530      if (usecprj==0) then
531        call pawcprj_free(cwaveprj_tmp)
532        ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
533      end if
534      nullify(cwaveprj_ptr)
535    else
536      cpopt=-1 ; paw_opt=0
537      call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
538 &     paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1_)
539    end if
540 
541 !  Electric field perturbation without Berry phase
542 !  -------------------------------------------
543  else if(ipert==natom+2 .and. &
544 &   (berryopt/=4 .and. berryopt/=6 .and. berryopt/=7 .and. &
545 &   berryopt/=14 .and. berryopt/=16 .and. berryopt/=17) .and.(optnl>0.or.sij_opt/=0))then
546 !  gvnl1 was already initialized in the calling routine, by reading a ddk file
547 !  It contains |i du^(0)/dk_band>
548 
549    if (gs_hamkq%usepaw==1) then
550      if (usecprj==1) then
551        cwaveprj_ptr => cwaveprj
552      else
553        ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor))
554        call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
555        cwaveprj_ptr => cwaveprj_tmp
556      end if
557      if (opt_gvnl1==2.and.optnl>=1) then
558 
559 !      PAW: Compute application of S^(0) to ddk WF
560        cpopt=-1 ; choice=1 ; paw_opt=3 ; signs=2
561        ABI_ALLOCATE(nonlop_out,(2,npw1*my_nspinor))
562        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,&
563 &       paw_opt,signs,nonlop_out,tim_nonlop,gvnl1_,vectout_dum)
564 !$OMP PARALLEL DO
565        do ipw=1,npw1*my_nspinor
566          gvnl1_(:,ipw)=nonlop_out(:,ipw)
567        end do
568 
569 !      PAW: Compute part of H^(1) due to derivative of S
570        cpopt=4*usecprj ; choice=51 ; paw_opt=3 ; signs=2
571        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,idir,(/lambda/),mpi_enreg,1,nnlout,&
572 &       paw_opt,signs,nonlop_out,tim_nonlop,cwave,vectout_dum)
573        if(compute_conjugate) then
574 !$OMP PARALLEL DO
575          do ipw=1,npw1*my_nspinor ! Note the multiplication by -i
576            gvnl1_(1,ipw)=gvnl1_(1,ipw)+nonlop_out(2,ipw)
577            gvnl1_(2,ipw)=gvnl1_(2,ipw)-nonlop_out(1,ipw)
578          end do
579        else
580 !$OMP PARALLEL DO
581          do ipw=1,npw1*my_nspinor ! Note the multiplication by i
582            gvnl1_(1,ipw)=gvnl1_(1,ipw)-nonlop_out(2,ipw)
583            gvnl1_(2,ipw)=gvnl1_(2,ipw)+nonlop_out(1,ipw)
584          end do
585        end if
586 
587 !      PAW: Compute part of H^(1) due to derivative of electric field part of Dij
588        cpopt=2 ; choice=1 ; paw_opt=1 ; signs=2
589        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,&
590 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr)
591 !$OMP PARALLEL DO
592        do ipw=1,npw1*my_nspinor
593          gvnl1_(:,ipw)=gvnl1_(:,ipw)+nonlop_out(:,ipw)
594        end do
595        ABI_DEALLOCATE(nonlop_out)
596 
597      end if ! opt_gvnl1==2
598 
599 !    PAW: Compute derivatives due to part of D_ij^(1) depending on VHxc^(1)
600      if (optnl>=2) then
601        ABI_ALLOCATE(gvnl2,(2,npw1*my_nspinor))
602        cpopt=-1+3*usecprj;if (opt_gvnl1==2) cpopt=2
603        choice=1 ; paw_opt=1 ; signs=2
604        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,0,(/lambda/),mpi_enreg,1,nnlout,&
605 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc)
606      end if
607 
608      if (sij_opt==1) then
609 !$OMP PARALLEL DO
610        do ipw=1,npw1*my_nspinor
611          gs1c(:,ipw)=zero
612        end do
613      end if
614      if (usecprj==0) then
615        call pawcprj_free(cwaveprj_tmp)
616        ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
617      end if
618      nullify(cwaveprj_ptr)
619    end if  ! PAW
620 
621 !  Electric field perturbation with Berry phase
622 !  -------------------------------------------
623  else if(ipert==natom+2 .and. &
624 &   (berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. &
625 &   berryopt==14 .or. berryopt==16 .or. berryopt==17 ) .and.(optnl>0.or.sij_opt/=0))then
626 
627    if (optnl>=1) then
628      do ipw=1,npw1*my_nspinor
629        gvnl1_(1,ipw)=-grad_berry(2,ipw)
630        gvnl1_(2,ipw)= grad_berry(1,ipw)
631      end do
632    end if
633    if (sij_opt==1) then
634 !$OMP PARALLEL DO
635      do ipw=1,npw1*my_nspinor
636        gs1c(:,ipw)=zero
637      end do
638    end if
639 
640 !  Strain perturbation
641 !  -------------------------------------------
642  else if ((ipert==natom+3.or.ipert==natom+4).and.(optnl>0.or.sij_opt/=0)) then
643 
644    istr=idir;if(ipert==natom+4) istr=istr+3
645 
646 !  PAW:
647    if (gs_hamkq%usepaw==1) then
648 
649      if (usecprj==1) then
650        cwaveprj_ptr => cwaveprj
651      else
652        ABI_DATATYPE_ALLOCATE(cwaveprj_tmp,(natom,my_nspinor))
653        call pawcprj_alloc(cwaveprj_tmp,1,gs_hamkq%dimcprj)
654        cwaveprj_ptr => cwaveprj_tmp
655      end if
656 
657 !    1- Compute derivatives due to projectors |p_i>^(1)
658 !    All atoms contribute
659      cpopt=-1+5*usecprj ; choice=3 ; signs=2
660      paw_opt=1;if (sij_opt/=0) paw_opt=sij_opt+3
661      call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,&
662 &     paw_opt,signs,gs1c,tim_nonlop,cwave,gvnl1_)
663 
664 !    2- Compute derivatives due to frozen part of D_ij^(1) (independent of VHxc^(1))
665 !    All atoms contribute
666      if (optnl>=1) then
667        ABI_ALLOCATE(nonlop_out,(2,npw1*my_nspinor))
668        cpopt=1+3*usecprj ; choice=1 ; signs=2 ; paw_opt=1
669        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,&
670 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave,nonlop_out,enl=rf_hamkq%e1kbfr)
671 !$OMP PARALLEL DO
672        do ipw=1,npw1*my_nspinor
673          gvnl1_(:,ipw)=gvnl1_(:,ipw)+nonlop_out(:,ipw)
674        end do
675        ABI_DEALLOCATE(nonlop_out)
676      end if
677 
678 !    3- Compute derivatives due to part of D_ij^(1) depending on VHxc^(1)
679 !    All atoms contribute
680      if (optnl>=2) then
681        ABI_ALLOCATE(gvnl2,(2,npw1*my_nspinor))
682        cpopt=4 ; choice=1 ; signs=2 ; paw_opt=1
683        call nonlop(choice,cpopt,cwaveprj_ptr,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,&
684 &       paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl2,enl=rf_hamkq%e1kbsc)
685      end if
686 
687      if (usecprj==0) then
688        call pawcprj_free(cwaveprj_tmp)
689        ABI_DATATYPE_DEALLOCATE(cwaveprj_tmp)
690      end if
691      nullify(cwaveprj_ptr)
692 
693 !    Norm-conserving psps:
694    else
695 !    Compute only derivatives due to projectors |p_i>^(1)
696      choice=3 ; cpopt=-1 ; signs=2 ; paw_opt=0
697      call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamkq,istr,(/lambda/),mpi_enreg,1,nnlout,&
698 &     paw_opt,signs,svectout_dum,tim_nonlop,cwave,gvnl1_)
699      if (sij_opt==1) then
700 !$OMP PARALLEL DO
701        do ipw=1,npw1*my_nspinor
702          gs1c(:,ipw)=zero
703        end do
704      end if
705    end if
706 
707 !  No non-local part
708 !  -------------------------------------------
709  else if (usevnl>0.or.(sij_opt/=0)) then
710 
711    if (optnl>=1) then
712 !$OMP PARALLEL DO
713      do ipw=1,npw1*my_nspinor
714        gvnl1_(:,ipw)=zero
715      end do
716    end if
717    if (sij_opt/=0) then
718 !$OMP PARALLEL DO
719      do ipw=1,npw1*my_nspinor
720        gs1c(:,ipw)=zero
721      end do
722    end if
723 
724  end if
725 
726 !======================================================================
727 !== Apply the 1st-order kinetic operator to the wavefunction
728 !== (add it to nl contribution)
729 !======================================================================
730 
731 !Phonon perturbation or Electric field perturbation
732 !-------------------------------------------
733 !No kinetic contribution
734 
735 !k-point perturbation or Strain perturbation
736 !-------------------------------------------
737 
738  usevnl2=allocated(gvnl2)
739  has_kin=(ipert==natom+1.or.ipert==natom+3.or.ipert==natom+4)
740  if (associated(gs_hamkq%kinpw_kp)) then
741    kinpw1 => gs_hamkq%kinpw_kp
742  else if (optnl>=1.or.usevnl2.or.has_kin) then
743    msg='need kinpw1 allocated!'
744    MSG_BUG(msg)
745  end if
746  if (associated(rf_hamkq%dkinpw_k)) then
747    dkinpw => rf_hamkq%dkinpw_k
748  else if (has_kin) then
749    msg='need dkinpw allocated!'
750    MSG_BUG(msg)
751  end if
752 
753  if (has_kin) then
754 !  Remember that npw=npw1 for ddk perturbation
755    do ispinor=1,my_nspinor
756 !$OMP PARALLEL DO PRIVATE(ipw,ipws) SHARED(cwave,ispinor,gvnl1_,dkinpw,kinpw1,npw,my_nspinor)
757      do ipw=1,npw
758        ipws=ipw+npw*(ispinor-1)
759        if(kinpw1(ipw)<huge(zero)*1.d-11)then
760          gvnl1_(1,ipws)=gvnl1_(1,ipws)+dkinpw(ipw)*cwave(1,ipws)
761          gvnl1_(2,ipws)=gvnl1_(2,ipws)+dkinpw(ipw)*cwave(2,ipws)
762        else
763          gvnl1_(1,ipws)=zero
764          gvnl1_(2,ipws)=zero
765        end if
766      end do
767    end do
768  end if
769 
770 !======================================================================
771 !== Sum contributions to get the application of H^(1) to the wf
772 !======================================================================
773 !Also filter the wavefunctions for large modified kinetic energy
774 
775 !Add non-local+kinetic to local part
776  if (optnl>=1.or.has_kin) then
777    do ispinor=1,my_nspinor
778      ipws=(ispinor-1)*npw1
779 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1c,gvnl1_,kinpw1,ipws,npw1)
780      do ipw=1+ipws,npw1+ipws
781        if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then
782          gh1c(1,ipw)=gh1c(1,ipw)+gvnl1_(1,ipw)
783          gh1c(2,ipw)=gh1c(2,ipw)+gvnl1_(2,ipw)
784        else
785          gh1c(1,ipw)=zero
786          gh1c(2,ipw)=zero
787        end if
788      end do
789    end do
790  end if
791 
792 !PAW: add non-local part due to first order change of VHxc
793  if (usevnl2) then
794    do ispinor=1,my_nspinor
795      ipws=(ispinor-1)*npw1
796 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gh1c,gvnl2,kinpw1,ipws,npw1)
797      do ipw=1+ipws,npw1+ipws
798        if(kinpw1(ipw-ipws)<huge(zero)*1.d-11)then
799          gh1c(1,ipw)=gh1c(1,ipw)+gvnl2(1,ipw)
800          gh1c(2,ipw)=gh1c(2,ipw)+gvnl2(2,ipw)
801        end if
802      end do
803    end do
804    ABI_DEALLOCATE(gvnl2)
805  end if
806 
807  if (usevnl==1) then
808    nullify(gvnl1_)
809  else
810    ABI_DEALLOCATE(gvnl1_)
811  end if
812 
813  call timab(196+tim_getgh1c,1,tsec)
814 
815  DBG_EXIT("COLL")
816 
817 end subroutine getgh1c

ABINIT/m_getgh1c [ Modules ]

[ Top ] [ Modules ]

NAME

  m_getgh1c

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (XG, DRH, MT, SPr)
  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_getgh1c
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32 
33  use defs_abitypes, only : MPI_type, dataset_type
34  use defs_datatypes, only : pseudopotential_type
35  use m_time,        only : timab
36  use m_pawcprj,     only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy, pawcprj_lincom, pawcprj_axpby
37  use m_kg,          only : kpgstr, mkkin, mkkpg
38  use m_mkffnl,      only : mkffnl
39  use m_pawfgr,      only : pawfgr_type
40  use m_fft,         only : fftpac, fourwf
41  use m_hamiltonian, only : gs_hamiltonian_type, rf_hamiltonian_type,&
42 &                          load_k_hamiltonian, load_kprime_hamiltonian,&
43 &                          load_k_rf_hamiltonian
44  use m_cgtools,          only : projbd
45  use m_nonlop,           only : nonlop
46  use m_fourier_interpol, only : transgrid
47 
48  implicit none
49 
50  private

m_hamiltonian/getgh1c_setup [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  getgh1c_setup

FUNCTION

INPUTS

OUTPUT

PARENTS

      dfpt_vtorho,m_gkk,m_phgamma,m_phpi,m_sigmaph

CHILDREN

      kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian,load_kprime_hamiltonian
      mkffnl,mkkin,mkkpg

SOURCE

 963 subroutine getgh1c_setup(gs_hamkq,rf_hamkq,dtset,psps,kpoint,kpq,idir,ipert,&           ! In
 964 &                natom,rmet,gprimd,gmet,istwf_k,npw_k,npw1_k,&                          ! In
 965 &                useylmgr1,kg_k,ylm_k,kg1_k,ylm1_k,ylmgr1_k,dkinpw,nkpg,&               ! In
 966 &                nkpg1,kpg_k,kpg1_k,kinpw1,ffnlk,ffnl1,ph3d,ph3d1,&                     ! Out
 967 &                ddkinpw,dkinpw2,rf_hamk_dir2,ffnl1_test)                               ! Optional
 968 
 969 
 970 !This section has been created automatically by the script Abilint (TD).
 971 !Do not modify the following lines by hand.
 972 #undef ABI_FUNC
 973 #define ABI_FUNC 'getgh1c_setup'
 974 !End of the abilint section
 975 
 976  implicit none
 977 
 978 !Arguments ------------------------------------
 979 !scalars
 980  integer,intent(in) :: idir,ipert,istwf_k,npw_k,npw1_k,natom,useylmgr1
 981  integer,intent(out) :: nkpg,nkpg1
 982  type(gs_hamiltonian_type),intent(inout) :: gs_hamkq
 983  type(rf_hamiltonian_type),intent(inout) :: rf_hamkq
 984  type(rf_hamiltonian_type),intent(inout),optional :: rf_hamk_dir2
 985  type(dataset_type),intent(in) :: dtset
 986  type(pseudopotential_type),intent(in) :: psps
 987 !arrays
 988  integer,intent(in) :: kg_k(3,npw_k),kg1_k(3,npw1_k)
 989  real(dp),intent(in) :: kpoint(3),kpq(3),gmet(3,3),gprimd(3,3),rmet(3,3)
 990  real(dp),intent(in) :: ylm_k(npw_k,psps%mpsang*psps%mpsang*psps%useylm)
 991  real(dp),intent(in) :: ylmgr1_k(npw1_k,3+6*((ipert-natom)/10),psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)
 992  real(dp),intent(in) :: ylm1_k(npw1_k,psps%mpsang*psps%mpsang*psps%useylm)
 993  real(dp),allocatable,intent(out) :: dkinpw(:),kinpw1(:),ffnlk(:,:,:,:),ffnl1(:,:,:,:)
 994  real(dp),allocatable,intent(out),optional :: dkinpw2(:),ddkinpw(:),ffnl1_test(:,:,:,:)
 995  real(dp),allocatable,intent(out) :: kpg_k(:,:),kpg1_k(:,:),ph3d(:,:,:),ph3d1(:,:,:)
 996 
 997 !Local variables-------------------------------
 998 !scalars
 999  integer :: dimffnl1,dimffnlk,ider,idir0,idir1,idir2,istr,ntypat,print_info
1000  logical :: qne0
1001 !arrays
1002  real(dp) :: ylmgr_dum(1,1,1)
1003 
1004 ! *************************************************************************
1005 
1006  if(.not.present(ddkinpw) .and. ipert==natom+10) then
1007    MSG_BUG("ddkinpw is not optional for ipert=natom+10.")
1008  end if
1009  if(.not.present(dkinpw2) .and. ipert==natom+10 .and. idir>3) then
1010    MSG_BUG("dkinpw2 is not optional for ipert=natom+10 and idir>3.")
1011  end if
1012  if(.not.present(rf_hamk_dir2) .and. ((ipert==natom+10 .and. idir>3) .or. ipert==natom+11)) then
1013    MSG_BUG("rf_hamk_dir2 is not optional for ipert=natom+10 (with idir>3) or ipert=natom+11.")
1014  end if
1015 
1016  ntypat = psps%ntypat
1017  qne0=((kpq(1)-kpoint(1))**2+(kpq(2)-kpoint(2))**2+(kpq(3)-kpoint(3))**2>=tol14)
1018 
1019 !Compute (k+G) vectors
1020  nkpg=0;if(ipert>=1.and.ipert<=natom) nkpg=3*dtset%nloalg(3)
1021  ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
1022  if (nkpg>0) then
1023    call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1024  end if
1025 
1026 !Compute (k+q+G) vectors
1027  nkpg1=0;if(ipert>=1.and.ipert<=natom) nkpg1=3*dtset%nloalg(3)
1028  ABI_ALLOCATE(kpg1_k,(npw1_k,nkpg1))
1029  if (nkpg1>0) then
1030    call mkkpg(kg1_k,kpg1_k,kpq(:),nkpg1,npw1_k)
1031  end if
1032 
1033 !===== Preparation of the non-local contributions
1034 
1035  dimffnlk=0;if (ipert<=natom) dimffnlk=1
1036  ABI_ALLOCATE(ffnlk,(npw_k,dimffnlk,psps%lmnmax,ntypat))
1037 
1038 !Compute nonlocal form factors ffnlk at (k+G)
1039 !(only for atomic displacement perturbation)
1040  if (ipert<=natom) then
1041    ider=0;idir0=0
1042    call mkffnl(psps%dimekb,dimffnlk,psps%ekb,ffnlk,psps%ffspl,&
1043 &   gmet,gprimd,ider,idir0,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,&
1044 &   psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,npw_k,ntypat,&
1045 &   psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm_k,ylmgr_dum)
1046  end if
1047 
1048 !Compute nonlocal form factors ffnl1 at (k+q+G)
1049  !-- Atomic displacement perturbation
1050  if (ipert<=natom) then
1051    ider=0;idir0=0
1052  !-- k-point perturbation (1st-derivative)
1053  else if (ipert==natom+1) then
1054    ider=1;idir0=idir
1055  !-- k-point perturbation (2nd-derivative)
1056  else if (ipert==natom+10.or.ipert==natom+11) then
1057    ider=2;idir0=4
1058  !-- Electric field perturbation
1059  else if (ipert==natom+2) then
1060    if (psps%usepaw==1) then
1061      ider=1;idir0=idir
1062    else
1063      ider=0;idir0=0
1064    end if
1065  !-- Strain perturbation
1066  else if (ipert==natom+3.or.ipert==natom+4) then
1067    if (ipert==natom+3) istr=idir
1068    if (ipert==natom+4) istr=idir+3
1069    ider=1;idir0=-istr
1070  !-- Magnetic field perturbation ( SPr, Zeeman )
1071  else if(ipert==natom+5)then
1072    ider=0;idir0=0
1073  end if
1074 
1075 !Compute nonlocal form factors ffnl1 at (k+q+G), for all atoms
1076  dimffnl1=1+ider
1077  if (ider==1.and.idir0==0) dimffnl1=2+2*psps%useylm
1078  if (ider==2.and.idir0==4) dimffnl1=3+7*psps%useylm
1079  ABI_ALLOCATE(ffnl1,(npw1_k,dimffnl1,psps%lmnmax,ntypat))
1080  call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1,psps%ffspl,gmet,gprimd,ider,idir0,&
1081 & psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
1082 & npw1_k,ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
1083 
1084 !Compute ffnl for nonlop with signs = 1
1085  print_info = 0
1086  if (dtset%prtvol==-19.or.dtset%prtvol==-20.or.dtset%prtvol==-21.or.dtset%nonlinear_info>=3) then
1087    print_info = 1
1088  end if
1089  if (present(ffnl1_test).and.print_info/=0.and.(ipert==natom+10.or.ipert==natom+11)) then
1090    ABI_ALLOCATE(ffnl1_test,(npw1_k,dimffnl1,psps%lmnmax,psps%ntypat))
1091    idir0 = 0 ! for nonlop with signs = 1
1092    call mkffnl(psps%dimekb,dimffnl1,psps%ekb,ffnl1_test,psps%ffspl,gs_hamkq%gmet,gs_hamkq%gprimd,ider,idir0,&
1093 &   psps%indlmn,kg1_k,kpg1_k,kpq,psps%lmnmax,psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg1,&
1094 &   npw1_k,psps%ntypat,psps%pspso,psps%qgrid_ff,rmet,psps%usepaw,psps%useylm,ylm1_k,ylmgr1_k)
1095  end if
1096 
1097 !===== Preparation of the kinetic contributions
1098 
1099 !Note that not all these arrays should be allocated in the general case when wtk_k vanishes
1100 
1101 !Compute (1/2) (2 Pi)**2 (k+q+G)**2:
1102  ABI_ALLOCATE(kinpw1,(npw1_k))
1103  kinpw1(:)=zero
1104  call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg1_k,kinpw1,kpq,npw1_k,0,0)
1105 
1106  ABI_ALLOCATE(dkinpw,(npw_k)) ! 1st derivative (1st direction)
1107  dkinpw(:)=zero
1108  if(ipert==natom+10 .and. idir>3) then
1109    ABI_ALLOCATE(dkinpw2,(npw_k)) ! 1st derivative (2nd directions)
1110    dkinpw2(:)=zero
1111  end if
1112  if(ipert==natom+10) then
1113    ABI_ALLOCATE(ddkinpw,(npw_k)) ! 2nd derivative
1114    ddkinpw(:)=zero
1115  end if
1116 
1117 !-- k-point perturbation (1st-derivative)
1118  if (ipert==natom+1) then
1119 !  Compute the derivative of the kinetic operator vs k
1120    call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative
1121  end if
1122 
1123 !-- k-point perturbation (2nd-derivative)
1124  if (ipert==natom+10.or.ipert==natom+11) then
1125 !  Compute the derivative of the kinetic operator vs k in kinpw, second and first orders
1126    if(ipert==natom+10 .and. idir<=3) then
1127      call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir,0) ! 1st derivative
1128      call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir,idir) ! 2nd derivative
1129    else
1130      select case(idir)
1131 !      Diagonal terms :
1132      case(1)
1133        idir1 = 1
1134        idir2 = 1
1135      case(2)
1136        idir1 = 2
1137        idir2 = 2
1138      case(3)
1139        idir1 = 3
1140        idir2 = 3
1141 !      Upper triangular terms :
1142      case(4)
1143        idir1 = 2
1144        idir2 = 3
1145      case(5)
1146        idir1 = 1
1147        idir2 = 3
1148      case(6)
1149        idir1 = 1
1150        idir2 = 2
1151 !      Lower triangular terms :
1152      case(7)
1153        idir1 = 3
1154        idir2 = 2
1155      case(8)
1156        idir1 = 3
1157        idir2 = 1
1158      case(9)
1159        idir1 = 2
1160        idir2 = 1
1161      end select
1162      call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw,kpoint,npw_k,idir1,0) !  1st derivative, idir1
1163      if(ipert==natom+10) then
1164        call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,dkinpw2,kpoint,npw_k,idir2,0) ! 1st derivative, idir2
1165        call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,ddkinpw,kpoint,npw_k,idir1,idir2) ! 2nd derivative
1166      end if
1167    end if
1168  end if
1169 
1170  !-- Strain perturbation
1171  if (ipert==natom+3.or.ipert==natom+4) then
1172    if (ipert==natom+3) istr=idir
1173    if (ipert==natom+4) istr=idir+3
1174 !  Compute the derivative of the kinetic operator vs strain
1175    call kpgstr(dkinpw,dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,gprimd,istr,kg_k,kpoint,npw_k)
1176  end if
1177 
1178 !===== Load the k/k+q dependent parts of the Hamiltonian
1179 
1180 !Load k-dependent part in the Hamiltonian datastructure
1181  ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamkq%matblk))
1182  call load_k_hamiltonian(gs_hamkq,kpt_k=kpoint,npw_k=npw_k,istwf_k=istwf_k,kg_k=kg_k,kpg_k=kpg_k,&
1183 & ph3d_k=ph3d,compute_ph3d=.true.,compute_gbound=.true.)
1184  if (size(ffnlk)>0) then
1185    call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnlk)
1186  else
1187    call load_k_hamiltonian(gs_hamkq,ffnl_k=ffnl1)
1188  end if
1189 
1190 !Load k+q-dependent part in the Hamiltonian datastructure
1191 !    Note: istwf_k is imposed to 1 for RF calculations (should use istwf_kq instead)
1192  call load_kprime_hamiltonian(gs_hamkq,kpt_kp=kpq,npw_kp=npw1_k,istwf_kp=istwf_k,&
1193 & kinpw_kp=kinpw1,kg_kp=kg1_k,kpg_kp=kpg1_k,ffnl_kp=ffnl1,&
1194 & compute_gbound=.true.)
1195  if (qne0) then
1196    ABI_ALLOCATE(ph3d1,(2,npw1_k,gs_hamkq%matblk))
1197    call load_kprime_hamiltonian(gs_hamkq,ph3d_kp=ph3d1,compute_ph3d=.true.)
1198  end if
1199 
1200 !Load k-dependent part in the 1st-order Hamiltonian datastructure
1201  call load_k_rf_hamiltonian(rf_hamkq,npw_k=npw_k,dkinpw_k=dkinpw)
1202  if (ipert==natom+10) then
1203    call load_k_rf_hamiltonian(rf_hamkq,ddkinpw_k=ddkinpw)
1204    if (idir>3) then
1205      call load_k_rf_hamiltonian(rf_hamk_dir2,dkinpw_k=dkinpw2,ddkinpw_k=ddkinpw)
1206    end if
1207  end if
1208 
1209 end subroutine getgh1c_setup

m_hamiltonian/rf_transgrid_and_pack [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  rf_transgrid_and_pack

FUNCTION

 Set up local potential vlocal1 with proper dimensioning, from vtrial1
 taking into account the spin. Same thing for vlocal from vtrial.

INPUTS

  isppol=Spin index.
  nspden=Number of density components
  usepaw=1 if PAW, 0 for NC.
  cplex=1 if DFPT potential is real, 2 for complex
  nfftf=Number of FFT points on the FINE grid treated by this processor
  nfft=Number of FFT points on the COARSE grid treated by this processor
  ngfft(18)=Info on the coarse grid.
  pawfgr <type(pawfgr_type)>=fine grid parameters and related data
  mpi_enreg=information about MPI parallelization
  vtrial(nfftf,nspden)=GS Vtrial(r) on the DENSE mesh
  vtrial1(cplex*nfftf,nspden)=INPUT RF Vtrial(r) on the DENSE mesh

OUTPUT

  vlocal(n4,n5,n6,nvloc)= GS local potential in real space, on the augmented coarse fft grid
  vlocal1(cplex*n4,n5,n6,nvloc)= RF local potential in real space, on the augmented coarse fft grid

PARENTS

      dfpt_vtorho,m_gkk,m_phgamma,m_phpi,m_sigmaph

CHILDREN

      kpgstr,load_k_hamiltonian,load_k_rf_hamiltonian,load_kprime_hamiltonian
      mkffnl,mkkin,mkkpg

SOURCE

856 subroutine rf_transgrid_and_pack(isppol,nspden,usepaw,cplex,nfftf,nfft,ngfft,nvloc,&
857 &                                pawfgr,mpi_enreg,vtrial,vtrial1,vlocal,vlocal1)
858 
859 
860 !This section has been created automatically by the script Abilint (TD).
861 !Do not modify the following lines by hand.
862 #undef ABI_FUNC
863 #define ABI_FUNC 'rf_transgrid_and_pack'
864 !End of the abilint section
865 
866  implicit none
867 
868 !Arguments ------------------------------------
869 !scalars
870  integer,intent(in) :: isppol,nspden,usepaw,cplex,nfftf,nfft,nvloc
871  type(pawfgr_type),intent(in) :: pawfgr
872  type(MPI_type),intent(in) :: mpi_enreg
873 !arrays
874  integer,intent(in) :: ngfft(18)
875  real(dp),intent(in),target :: vtrial(nfftf,nspden)
876  real(dp),intent(inout),target :: vtrial1(cplex*nfftf,nspden)
877  real(dp),intent(out) :: vlocal(ngfft(4),ngfft(5),ngfft(6),nvloc)
878  real(dp),intent(out) :: vlocal1(cplex*ngfft(4),ngfft(5),ngfft(6),nvloc)
879 
880 !Local variables-------------------------------
881 !scalars
882  integer :: n1,n2,n3,n4,n5,n6,paral_kgb,ispden
883 !arrays
884  real(dp) :: rhodum(1)
885  real(dp), ABI_CONTIGUOUS pointer :: vtrial_ptr(:,:),vtrial1_ptr(:,:)
886  real(dp),allocatable :: cgrvtrial(:,:),cgrvtrial1(:,:),vlocal_tmp(:,:,:),vlocal1_tmp(:,:,:)
887 
888 ! *************************************************************************
889 
890  n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
891  n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
892  paral_kgb = mpi_enreg%paral_kgb
893 
894  if (nspden/=4) then
895    vtrial_ptr=>vtrial
896    if (usepaw==0.or.pawfgr%usefinegrid==0) then
897      call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal(:,:,:,1),2)
898      call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1,vlocal1(:,:,:,1),2)
899    else
900      ABI_ALLOCATE(cgrvtrial,(nfft,nspden))
901      call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr)
902      call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,cgrvtrial,vlocal(:,:,:,1),2)
903      ABI_DEALLOCATE(cgrvtrial)
904      ABI_ALLOCATE(cgrvtrial,(cplex*nfft,nspden))
905      call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial1)
906      call fftpac(isppol,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,cgrvtrial,vlocal1(:,:,:,1),2)
907      ABI_DEALLOCATE(cgrvtrial)
908    end if
909    nullify(vtrial_ptr)
910  else  ! nspden==4 non-collinear magnetism
911    vtrial_ptr=>vtrial
912    vtrial1_ptr=>vtrial1
913    ABI_ALLOCATE(vlocal_tmp,(n4,n5,n6))
914    ABI_ALLOCATE(vlocal1_tmp,(cplex*n4,n5,n6))
915    if (usepaw==0.or.pawfgr%usefinegrid==0) then
916      do ispden=1,nspden
917        call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2)
918        vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
919        call fftpac(ispden,mpi_enreg,nspden,cplex*n1,n2,n3,cplex*n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2)
920        vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:)
921      end do
922    else ! TODO FR EB check the correctness of the following lines for PAW calculations
923      ABI_ALLOCATE(cgrvtrial,(nfft,nspden))
924      ABI_ALLOCATE(cgrvtrial1,(nfft,nspden))
925      call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial_ptr)
926      call transgrid(cplex,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial1,vtrial1_ptr)
927      do ispden=1,nspden
928        call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial_ptr,vlocal_tmp,2)
929        vlocal(:,:,:,ispden)=vlocal_tmp(:,:,:)
930        call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial1_ptr,vlocal1_tmp,2)
931        vlocal1(:,:,:,ispden)=vlocal1_tmp(:,:,:)
932      end do
933      ABI_DEALLOCATE(cgrvtrial)
934    end if
935    ABI_DEALLOCATE(vlocal_tmp)
936    ABI_DEALLOCATE(vlocal1_tmp)
937  end if !nspden
938 
939 end subroutine rf_transgrid_and_pack