TABLE OF CONTENTS
ABINIT/dotprodm_vn [ Functions ]
NAME
dotprodm_vn
FUNCTION
For a set of densities and a set of potentials, compute the dot product (integral over FFT grid) of each pair, to obtain a series of 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 spin-polarisation (if nspden=2), or the magnetization vector (if nspden=4). Need the index of the first density/potential pair to be treated, in each array, and the number of pairs to be treated. Might be used to compute just one dot product, in a big array, such as to avoid copying the density and potential from a big array to a temporary place.
COPYRIGHT
Copyright (C) 1999-2018 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 cpldot=if 1, the dot array is real, if 2, the dot array is complex (not coded yet for nspden=4) denarr(cplex*nfft,nspden,nden)=real space density on FFT grid id=index of the first density to be treated in the denarr array ip=index of the first potential to be treated in the potarr array mpicomm=the mpi communicator used for the summation mpi_summarize=set it to .true. if parallelisation is done over FFT multd=number of densities to be treated multp=number of potentials to be treated nden=third dimension of the denarr array nfft= (effective) number of FFT grid points (for this processor) nfftot= total number of FFT grid points npot=third dimension of the potarr array nspden=number of spin-density components potarr(cplex*nfft,nspden,npot)=real space potential on FFT grid (will be complex conjugated if cplex=2 and cpldot=2) ucvol=unit cell volume (Bohr**3)
OUTPUT
dot(cpldot,multp,multd)= series of values of the dot product potential/density
SIDE EFFECTS
NOTES
Concerning storage when nspden=4:
cplex=1:
V are stored as : V^11, V^22, Re[V^12], Im[V^12] (complex, hermitian)
N are stored as : n, m_x, m_y, m_z (real)
cplex=2:
V are stored as : V^11, V^22, V^12, i.V^21 (complex)
N are stored as : n, m_x, m_y, mZ (complex)
PARENTS
aprxdr
CHILDREN
timab,xmpi_sum
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 78 subroutine dotprodm_vn(cplex,cpldot,denarr,dot,id,ip,mpicomm, mpi_summarize,multd,multp,& 79 & nden,nfft,nfftot,npot,nspden,potarr,ucvol) 80 81 use defs_basis 82 use m_errors 83 use m_profiling_abi 84 use m_xmpi 85 86 !This section has been created automatically by the script Abilint (TD). 87 !Do not modify the following lines by hand. 88 #undef ABI_FUNC 89 #define ABI_FUNC 'dotprodm_vn' 90 use interfaces_18_timing 91 !End of the abilint section 92 93 implicit none 94 95 !Arguments ------------------------------------ 96 !scalars 97 integer,intent(in) :: cpldot,cplex,id,ip,multd,multp,nden,nfft,nfftot,npot 98 integer,intent(in) :: nspden,mpicomm 99 logical, intent(in) :: mpi_summarize 100 real(dp),intent(in) :: ucvol 101 !arrays 102 real(dp),intent(in) :: denarr(cplex*nfft,nspden,nden) 103 real(dp),intent(in) :: potarr(cplex*nfft,nspden,npot) 104 real(dp),intent(out) :: dot(cpldot,multp,multd) 105 106 !Local variables------------------------------- 107 !scalars 108 integer :: i1,i2,ierr,ir,jr 109 real(dp) :: ai,ar,dim11,dim12,dim21,dim22,dim_dn,dim_up,dre11,dre12,dre21 110 real(dp) :: dre22,dre_dn,dre_up,factor,pim11,pim12,pim21,pim22,pim_dn,pim_up 111 real(dp) :: pre11,pre12,pre21,pre22,pre_dn,pre_up 112 !arrays 113 real(dp) :: tsec(2) 114 115 ! ************************************************************************* 116 117 !Real or complex inputs are coded 118 DBG_CHECK(ANY(cplex==(/1,2/)),"Wrong cplex") 119 120 !Real or complex outputs are coded 121 DBG_CHECK(ANY(cpldot==(/1,2/)),"Wrong cpldot") 122 DBG_CHECK(ANY(nspden==(/1,2,4/)),"Wrong nspden") 123 124 DBG_CHECK(id >= 1,'Wrong id') 125 DBG_CHECK(ip >= 1,'Wrong id') 126 127 DBG_CHECK(multd >= 1,"wrong multd") 128 DBG_CHECK(multp >= 1,"wrong multp") 129 130 DBG_CHECK(nden-id-multd >=-1,'nden-id-multd') 131 DBG_CHECK(npot-ip-multp >=-1,'npot-ip-multp') 132 133 if(nspden==1)then 134 135 if(cpldot==1 .or. cplex==1 )then 136 137 do i2=1,multd 138 do i1=1,multp 139 ar=zero 140 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar) 141 do ir=1,cplex*nfft 142 ar=ar + potarr(ir,1,ip+i1-1)*denarr(ir,1,id+i2-1) 143 end do 144 dot(1,i1,i2)=ar 145 end do ! i1 146 end do ! i2 147 148 else ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot 149 150 do i2=1,multd 151 do i1=1,multp 152 ar=zero ; ai=zero 153 !$OMP PARALLEL DO PRIVATE(ir,jr) SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai) 154 do ir=1,nfft 155 jr=2*ir 156 ar=ar + potarr(jr-1,1,ip+i1-1)*denarr(jr-1,1,id+i2-1) & 157 & + potarr(jr ,1,ip+i1-1)*denarr(jr ,1,id+i2-1) 158 ai=ai + potarr(jr-1,1,ip+i1-1)*denarr(jr ,1,id+i2-1) & 159 & - potarr(jr ,1,ip+i1-1)*denarr(jr-1,1,id+i2-1) 160 end do 161 dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai 162 end do ! i1 163 end do ! i2 164 165 end if 166 167 else if(nspden==2)then 168 169 if(cpldot==1 .or. cplex==1 )then 170 171 do i2=1,multd 172 do i1=1,multp 173 ar=zero 174 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar) 175 do ir=1,cplex*nfft 176 ar=ar + potarr(ir,1,ip+i1-1)* denarr(ir,2,id+i2-1) & ! This is the spin up contribution 177 & + potarr(ir,2,ip+i1-1)*(denarr(ir,1,id+i2-1)-denarr(ir,2,id+i2-1)) ! This is the spin down contribution 178 end do 179 dot(1,i1,i2)=ar 180 end do ! i1 181 end do ! i2 182 183 else ! cpldot==2 and cplex==2 : one builds the imaginary part, from complex den/pot 184 185 do i2=1,multd 186 do i1=1,multp 187 ar=zero ; ai=zero 188 !$OMP PARALLEL DO PRIVATE(ir,jr,dre_up,dim_up,dre_dn,dim_dn,pre_up,pim_up,pre_dn,pim_dn) & 189 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai) 190 do ir=1,nfft 191 jr=2*ir 192 193 dre_up=denarr(jr-1,2,id+i2-1) 194 dim_up=denarr(jr ,2,id+i2-1) 195 dre_dn=denarr(jr-1,1,id+i2-1)-dre_up 196 dim_dn=denarr(jr ,1,id+i2-1)-dim_up 197 198 pre_up=potarr(jr-1,1,ip+i1-1) 199 pim_up=potarr(jr ,1,ip+i1-1) 200 pre_dn=potarr(jr-1,2,ip+i1-1) 201 pim_dn=potarr(jr ,2,ip+i1-1) 202 203 ar=ar + pre_up * dre_up & 204 & + pim_up * dim_up & 205 & + pre_dn * dre_dn & 206 & + pim_dn * dim_dn 207 ai=ai + pre_up * dim_up & 208 & - pim_up * dre_up & 209 & + pre_dn * dim_dn & 210 & - pim_dn * dre_dn 211 212 end do 213 dot(1,i1,i2)=ar ; dot(2,i1,i2)=ai 214 end do ! i1 215 end do ! i2 216 217 end if 218 219 else if(nspden==4)then 220 ! \rho{\alpha,\beta} V^{\alpha,\beta} = 221 ! rho*(V^{11}+V^{22})/2$ 222 ! + m_x Re(V^{12})- m_y Im{V^{12}}+ m_z(V^{11}-V^{22})/2 223 if (cplex==1) then 224 do i2=1,multd 225 do i1=1,multp 226 ar=zero 227 !$OMP PARALLEL DO PRIVATE(ir) SHARED(id,i1,i2,ip,cplex,nfft,denarr,potarr) REDUCTION(+:ar) 228 do ir=1,cplex*nfft 229 ar=ar+(potarr(ir,1,ip+i1-1)+potarr(ir,2,ip+i1-1))*half*denarr(ir,1,id+i2-1)& ! This is the density contrib 230 & + potarr(ir,3,ip+i1-1) *denarr(ir,2,id+i2-1)& ! This is the m_x contrib 231 & - potarr(ir,4,ip+i1-1) *denarr(ir,3,id+i2-1)& ! This is the m_y contrib 232 & +(potarr(ir,1,ip+i1-1)-potarr(ir,2,ip+i1-1))*half*denarr(ir,4,id+i2-1) ! This is the m_z contrib 233 end do 234 dot(1,i1,i2)=ar 235 end do ! i1 236 end do ! i2 237 else ! cplex=2 238 ! Note concerning storage when cplex=2: 239 ! V are stored as : v^11, v^22, V^12, i.V^21 (each are complex) 240 ! N are stored as : n, m_x, m_y, mZ (each are complex) 241 if (cpldot==1) then 242 do i2=1,multd 243 do i1=1,multp 244 ar=zero ; ai=zero 245 !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre22,pim22,pre12,pim12) & 246 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar) 247 do ir=1,nfft 248 jr=2*ir 249 dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2)) 250 dim11=half*(denarr(jr ,1,id+i2)+denarr(jr-1,4,id+i2)) 251 dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2)) 252 dim22=half*(denarr(jr ,1,id+i2)-denarr(jr-1,4,id+i2)) 253 dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr ,3,id+i2)) 254 dim12=half*(denarr(jr ,2,id+i2)-denarr(jr-1,3,id+i2)) 255 dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr ,3,id+i2)) 256 dim21=half*(denarr(jr ,2,id+i2)+denarr(jr-1,3,id+i2)) 257 pre11= potarr(jr-1,1,ip+i1) 258 pim11= potarr(jr ,1,ip+i1) 259 pre22= potarr(jr-1,2,ip+i1) 260 pim22= potarr(jr ,2,ip+i1) 261 pre12= potarr(jr-1,3,ip+i1) 262 pim12= potarr(jr ,3,ip+i1) 263 pre21= potarr(jr ,4,ip+i1) 264 pim21=-potarr(jr-1,4,ip+i1) 265 ar=ar + pre11 * dre11 & 266 & + pim11 * dim11 & 267 & + pre22 * dre22 & 268 & + pim22 * dim22 & 269 & + pre12 * dre12 & 270 & + pim12 * dim12 & 271 & + pre21 * dre21 & 272 & + pim21 * dim21 273 end do 274 dot(1,i1,i2)=ar 275 end do ! i1 276 end do ! i2 277 else !cpldot=2 278 do i2=1,multd 279 do i1=1,multp 280 ar=zero ; ai=zero 281 !$OMP PARALLEL DO PRIVATE(ir,jr,dre11,dim11,dre22,dim22,dre12,dim12,pre11,pim11,pre12,pim12,pre22,pim22) & 282 !$OMP&SHARED(id,i1,i2,ip,nfft,denarr,potarr) REDUCTION(+:ar,ai) 283 do ir=1,nfft 284 jr=2*ir 285 dre11=half*(denarr(jr-1,1,id+i2)+denarr(jr-1,4,id+i2)) 286 dim11=half*(denarr(jr ,1,id+i2)+denarr(jr-1,4,id+i2)) 287 dre22=half*(denarr(jr-1,1,id+i2)-denarr(jr-1,4,id+i2)) 288 dim22=half*(denarr(jr ,1,id+i2)-denarr(jr-1,4,id+i2)) 289 dre12=half*(denarr(jr-1,2,id+i2)+denarr(jr ,3,id+i2)) 290 dim12=half*(denarr(jr ,2,id+i2)-denarr(jr-1,3,id+i2)) 291 dre21=half*(denarr(jr-1,2,id+i2)-denarr(jr ,3,id+i2)) 292 dim21=half*(denarr(jr ,2,id+i2)+denarr(jr-1,3,id+i2)) 293 pre11= potarr(jr-1,1,ip+i1) 294 pim11= potarr(jr ,1,ip+i1) 295 pre22= potarr(jr-1,2,ip+i1) 296 pim22= potarr(jr ,2,ip+i1) 297 pre12= potarr(jr-1,3,ip+i1) 298 pim12= potarr(jr ,3,ip+i1) 299 pre21= potarr(jr ,4,ip+i1) 300 pim21=-potarr(jr-1,4,ip+i1) 301 ar=ar + pre11 * dre11 & 302 & + pim11 * dim11 & 303 & + pre22 * dre22 & 304 & + pim22 * dim22 & 305 & + pre12 * dre12 & 306 & + pim12 * dim12 & 307 & + pre21 * dre21 & 308 & + pim21 * dim21 309 ai=ai + pre11 * dim11 & 310 & - pim11 * dre11 & 311 & + pre22 * dim22 & 312 & - pim22 * dre22 & 313 & + pre12 * dim12 & 314 & - pim12 * dre12 & 315 & + pre21 * dim21 & 316 & - pim21 * dre21 317 end do 318 dot(1,i1,i2)=ar 319 dot(2,i1,i2)=ai 320 end do ! i1 321 end do ! i2 322 end if ! cpldot 323 end if ! cplex 324 end if ! nspden 325 326 factor=ucvol/dble(nfftot) 327 dot(:,:,:)=factor*dot(:,:,:) 328 329 !XG030513 : MPIWF reduction (addition) on dot is needed here 330 if (mpi_summarize) then 331 call timab(48,1,tsec) 332 call xmpi_sum(dot,mpicomm ,ierr) 333 call timab(48,2,tsec) 334 end if 335 336 if(cpldot==2 .and. cplex==1)dot(2,:,:)=zero 337 338 end subroutine dotprodm_vn