TABLE OF CONTENTS
ABINIT/dotprod_vn [ Functions ]
NAME
dotprod_vn
FUNCTION
Compute dot product of potential and density (integral over FFT grid), to obtain an energy-like quantity (so the usual dotproduct is divided by the number of FFT points, and multiplied by the primitive cell volume). Take into account the spin components of the density and potentials (nspden), and sum correctly over them. Note that the storage of densities and potentials is different : for potential, one stores the matrix components, while for the density, one stores the trace, and then, either the up component (if nspden=2), or the magnetization vector (if nspden=4).
COPYRIGHT
Copyright (C) 1999-2017 ABINIT group (XG) 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
cplex=if 1, real space functions on FFT grid are REAL, if 2, COMPLEX dens(cplex*nfft,nspden)=real space density on FFT grid mpi_enreg=informations about MPI parallelization nfft= (effective) number of FFT grid points (for this processor) nfftot= total number of FFT grid points nspden=number of spin-density components option= if 1, only the real part is computed if 2, both real and imaginary parts are computed (not yet coded) pot(cplex*nfft,nspden)=real space potential on FFT grid (will be complex conjugated if cplex=2 and option=2) ucvol=unit cell volume (Bohr**3)
OUTPUT
doti= imaginary part of the dot product, output only if option=2 (and cplex=2). dotr= real part
SIDE EFFECTS
NOTES
Concerning storage when nspden=4:
cplex=1:
V is stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
N is stored as : n, m_x, m_y, m_z (real)
cplex=2:
V is stored as : V^11, V^22, V^12, i.V^21 (complex)
N is stored as : n, m_x, m_y, mZ (complex)
PARENTS
dfpt_dyxc1,dfpt_eltfrxc,dfpt_nselt,dfpt_nstdy,dfpt_nstpaw,dfpt_rhotov dfptnl_loop,energy,newfermie1,odamix,prcrskerker2,rhohxc,rhotov,setvtr
CHILDREN
xmpi_sum
SOURCE
61 #if defined HAVE_CONFIG_H 62 #include "config.h" 63 #endif 64 65 66 #include "abi_common.h" 67 68 subroutine dotprod_vn(cplex,dens,dotr,doti,nfft,nfftot,nspden,option,pot,& 69 & ucvol,mpi_comm_sphgrid) 70 71 use defs_basis 72 use m_profiling_abi 73 use m_errors 74 75 use m_xmpi, only: xmpi_comm_size,xmpi_sum 76 77 !This section has been created automatically by the script Abilint (TD). 78 !Do not modify the following lines by hand. 79 #undef ABI_FUNC 80 #define ABI_FUNC 'dotprod_vn' 81 !End of the abilint section 82 83 implicit none 84 85 !Arguments ------------------------------------ 86 !scalars 87 integer,intent(in) :: cplex,nfft,nfftot,nspden,option 88 integer,intent(in),optional :: mpi_comm_sphgrid 89 real(dp),intent(in) :: ucvol 90 real(dp),intent(out) :: doti,dotr 91 !arrays 92 real(dp),intent(in) :: dens(cplex*nfft,nspden),pot(cplex*nfft,nspden) 93 94 !Local variables------------------------------- 95 !scalars 96 integer :: ierr,ifft,jfft 97 real(dp) :: dim11,dim12,dim21,dim22,dim_dn,dim_up,dre11,dre12,dre21,dre22 98 real(dp) :: dre_dn,dre_up,factor,nproc_sphgrid,pim11,pim12,pim21,pim22,pim_dn,pim_up,pre11 99 real(dp) :: pre12,pre21,pre22,pre_dn,pre_up 100 !arrays 101 real(dp) :: buffer2(2) 102 103 ! ************************************************************************* 104 105 !Real or complex inputs are coded 106 DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex") 107 DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden") 108 109 !Real or complex output are coded 110 DBG_CHECK(ANY(option==(/1,2/)),"Wrong option") 111 112 dotr=zero; doti=zero 113 114 if(nspden==1)then 115 116 if(option==1 .or. cplex==1 )then 117 !$OMP PARALLEL DO REDUCTION(+:dotr) 118 do ifft=1,cplex*nfft 119 dotr=dotr + pot(ifft,1)*dens(ifft,1) 120 end do 121 ! dotr = ddot(cplex*nfft,pot,1,dens,1) 122 123 else ! option==2 and cplex==2 : one builds the imaginary part, from complex den/pot 124 125 !$OMP PARALLEL DO PRIVATE(jfft) REDUCTION(+:dotr,doti) 126 do ifft=1,nfft 127 jfft=2*ifft 128 dotr=dotr + pot(jfft-1,1)*dens(jfft-1,1) + pot(jfft,1)*dens(jfft ,1) 129 doti=doti + pot(jfft-1,1)*dens(jfft ,1) - pot(jfft,1)*dens(jfft-1,1) 130 end do 131 132 end if 133 134 else if(nspden==2)then 135 136 if(option==1 .or. cplex==1 )then 137 !$OMP PARALLEL DO REDUCTION(+:dotr) 138 do ifft=1,cplex*nfft 139 dotr=dotr + pot(ifft,1)* dens(ifft,2) & ! This is the spin up contribution 140 & + pot(ifft,2)*(dens(ifft,1)-dens(ifft,2)) ! This is the spin down contribution 141 end do 142 143 else ! option==2 and cplex==2 : one builds the imaginary part, from complex den/pot 144 145 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr,doti) 146 do ifft=1,nfft 147 148 jfft=2*ifft 149 dre_up=dens(jfft-1,2) 150 dim_up=dens(jfft ,2) 151 dre_dn=dens(jfft-1,1)-dre_up 152 dim_dn=dens(jfft ,1)-dim_up 153 pre_up=pot(jfft-1,1) 154 pim_up=pot(jfft ,1) 155 pre_dn=pot(jfft-1,2) 156 pim_dn=pot(jfft ,2) 157 158 dotr=dotr + pre_up * dre_up & 159 & + pim_up * dim_up & 160 & + pre_dn * dre_dn & 161 & + pim_dn * dim_dn 162 doti=doti + pre_up * dim_up & 163 & - pim_up * dre_up & 164 & + pre_dn * dim_dn & 165 & - pim_dn * dre_dn 166 167 end do 168 end if 169 170 else if(nspden==4)then 171 ! \rho{\alpha,\beta} V^{\alpha,\beta} = 172 ! rho*(V^{11}+V^{22})/2$ 173 ! + m_x Re(V^{12})- m_y Im{V^{12}}+ m_z(V^{11}-V^{22})/2 174 if (cplex==1) then 175 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(nfft,dens,pot) REDUCTION(+:dotr) 176 do ifft=1,nfft 177 dotr=dotr + & 178 & (pot(ifft,1) + pot(ifft,2))*half*dens(ifft,1) & ! This is the density contrib 179 & + pot(ifft,3) *dens(ifft,2) & ! This is the m_x contrib 180 & - pot(ifft,4) *dens(ifft,3) & ! This is the m_y contrib 181 & +(pot(ifft,1) - pot(ifft,2))*half*dens(ifft,4) ! This is the m_z contrib 182 end do 183 else ! cplex=2 184 ! Note concerning storage when cplex=2: 185 ! V is stored as : v^11, v^22, V^12, i.V^21 (each are complex) 186 ! N is stored as : n, m_x, m_y, mZ (each are complex) 187 if (option==1) then 188 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr) 189 do ifft=1,nfft 190 jfft=2*ifft 191 dre11=half*(dens(jfft-1,1)+dens(jfft-1,4)) 192 dim11=half*(dens(jfft ,1)+dens(jfft-1,4)) 193 dre22=half*(dens(jfft-1,1)-dens(jfft-1,4)) 194 dim22=half*(dens(jfft ,1)-dens(jfft-1,4)) 195 dre12=half*(dens(jfft-1,2)+dens(jfft ,3)) 196 dim12=half*(dens(jfft ,2)-dens(jfft-1,3)) 197 dre21=half*(dens(jfft-1,2)-dens(jfft ,3)) 198 dim21=half*(dens(jfft ,2)+dens(jfft-1,3)) 199 pre11= pot(jfft-1,1) 200 pim11= pot(jfft ,1) 201 pre22= pot(jfft-1,2) 202 pim22= pot(jfft ,2) 203 pre12= pot(jfft-1,3) 204 pim12= pot(jfft ,3) 205 pre21= pot(jfft ,4) 206 pim21=-pot(jfft-1,4) 207 dotr=dotr + pre11 * dre11 & 208 & + pim11 * dim11 & 209 & + pre22 * dre22 & 210 & + pim22 * dim22 & 211 & + pre12 * dre12 & 212 & + pim12 * dim12 & 213 & + pre21 * dre21 & 214 & + pim21 * dim21 215 end do 216 else ! option=2 217 !$OMP PARALLEL DO DEFAULT(PRIVATE) SHARED(nfft,dens,pot) REDUCTION(+:dotr,doti) 218 do ifft=1,nfft 219 jfft=2*ifft 220 dre11=half*(dens(jfft-1,1)+dens(jfft-1,4)) 221 dim11=half*(dens(jfft ,1)+dens(jfft-1,4)) 222 dre22=half*(dens(jfft-1,1)-dens(jfft-1,4)) 223 dim22=half*(dens(jfft ,1)-dens(jfft-1,4)) 224 dre12=half*(dens(jfft-1,2)+dens(jfft ,3)) 225 dim12=half*(dens(jfft ,2)-dens(jfft-1,3)) 226 dre21=half*(dens(jfft-1,2)-dens(jfft ,3)) 227 dim21=half*(dens(jfft ,2)+dens(jfft-1,3)) 228 pre11= pot(jfft-1,1) 229 pim11= pot(jfft ,1) 230 pre22= pot(jfft-1,2) 231 pim22= pot(jfft ,2) 232 pre12= pot(jfft-1,3) 233 pim12= pot(jfft ,3) 234 pre21= pot(jfft ,4) 235 pim21=-pot(jfft-1,4) 236 dotr=dotr + pre11 * dre11 & 237 & + pim11 * dim11 & 238 & + pre22 * dre22 & 239 & + pim22 * dim22 & 240 & + pre12 * dre12 & 241 & + pim12 * dim12 & 242 & + pre21 * dre21 & 243 & + pim21 * dim21 244 doti=doti + pre11 * dim11 & 245 & - pim11 * dre11 & 246 & + pre22 * dim22 & 247 & - pim22 * dre22 & 248 & + pre12 * dim12 & 249 & - pim12 * dre12 & 250 & + pre21 * dim21 & 251 & - pim21 * dre21 252 end do 253 end if ! option 254 end if ! cplex 255 end if ! nspden 256 257 factor=ucvol/dble(nfftot) 258 dotr=factor*dotr 259 doti=factor*doti 260 261 !MPIWF reduction (addition) on dotr, doti is needed here 262 if(present(mpi_comm_sphgrid)) then 263 nproc_sphgrid=xmpi_comm_size(mpi_comm_sphgrid) 264 if(nproc_sphgrid>1) then 265 buffer2(1)=dotr 266 buffer2(2)=doti 267 call xmpi_sum(buffer2,mpi_comm_sphgrid,ierr) 268 dotr=buffer2(1) 269 doti=buffer2(2) 270 end if 271 end if 272 273 end subroutine dotprod_vn