1!! gen1int: compute one-electron integrals using rotational London atomic-orbitals 2!! Copyright 2009-2012 Bin Gao, and Andreas Thorvaldsen 3!! 4!! gen1int is free software: you can redistribute it and/or modify 5!! it under the terms of the GNU Lesser General Public License as published by 6!! the Free Software Foundation, either version 3 of the License, or 7!! (at your option) any later version. 8!! 9!! gen1int is distributed in the hope that it will be useful, 10!! but WITHOUT ANY WARRANTY; without even the implied warranty of 11!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12!! GNU Lesser General Public License for more details. 13!! 14!! You should have received a copy of the GNU Lesser General Public License 15!! along with gen1int. If not, see <http://www.gnu.org/licenses/>. 16!! 17!! This file contains the recurrence relations of nuclear attraction potential 18!! (multiplied by Cartesian multipole moment) integrals using primitive Hermite 19!! Gaussians. 20!! 21!! 2011-07-27, Bin Gao 22!! * first version 23 24#include "stdout.h" 25 26 !> \brief recurrence relations of nuclear attraction potential (multiplied by 27 !> Cartesian multipole moment) integrals using primitive Hermite Gaussians 28 !> \author Bin Gao 29 !> \date 2011-07-27 30 !> \param orders_hgto_bra is the range of orders of Hermite Gaussians on bra center 31 !> \param coord_bra contains the coordinates of bra center 32 !> \param exponent_bra is the exponent of primitive Gaussian of bra center 33 !> \param orders_hgto_ket is the range of orders of Hermite Gaussians on ket center 34 !> \param coord_ket contains the coordinates of ket center 35 !> \param exponent_ket is the exponent of primitive Gaussian of ket center 36 !> \param nucpot_origin contains the coordinates of nuclear potential origin 37 !> \param dipole_origin contains the coordinates of dipole origin 38 !> \param scal_const is the scale constant for nuclear attraction potential 39 !> \param order_geo_pot is the order of geometric derivatives on potential origin 40 !> \param order_mom is the order of Cartesian multipole moments 41 !> \param order_elec is the order of electronic derivatives 42 !> \param dim_hgto_bra is the dimension of Hermite Gaussians on bra center 43 !> \param dim_hgto_ket is the dimension of Hermite Gaussians on ket center 44 !> \param num_elec is the number of xyz components of electronic derivatives, 45 !> equals to \f$(\var(order_elec)+1)(\var(order_elec)+2)/2\f$ 46 !> \param num_mom is the number of xyz components of Cartesian multipole moment, 47 !> equals to \f$(\var(order_mom)+1)(\var(order_mom)+2)/2\f$ 48 !> \param num_geo_pot is the number of geometric derivatives on potential origin, 49 !> equals to \f$(\var(order_geo_pot)+1)(\var(order_geo_pot)+2)/2\f$ 50 !> \return hgto_pints contains the primitive HGTO integrals 51 subroutine prim_hgto_nucpot(orders_hgto_bra, coord_bra, exponent_bra, & 52 orders_hgto_ket, coord_ket, exponent_ket, & 53 order_elec, nucpot_origin, dipole_origin, & 54 scal_const, order_mom, order_geo_pot, & 55 dim_hgto_bra, dim_hgto_ket, num_elec, & 56 num_mom, num_geo_pot, hgto_pints) 57 use xkind 58 implicit none 59 integer, intent(in) :: orders_hgto_bra(2) 60 real(REALK), intent(in) :: coord_bra(3) 61 real(REALK), intent(in) :: exponent_bra 62 integer, intent(in) :: orders_hgto_ket(2) 63 real(REALK), intent(in) :: coord_ket(3) 64 real(REALK), intent(in) :: exponent_ket 65 integer, intent(in) :: order_elec 66 real(REALK), intent(in) :: nucpot_origin(3) 67 real(REALK), intent(in) :: dipole_origin(3) 68 real(REALK), intent(in) :: scal_const 69 integer, intent(in) :: order_mom 70 integer, intent(in) :: order_geo_pot 71 integer, intent(in) :: dim_hgto_bra 72 integer, intent(in) :: dim_hgto_ket 73 integer, intent(in) :: num_elec 74 integer, intent(in) :: num_mom 75 integer, intent(in) :: num_geo_pot 76 real(REALK), intent(out) :: hgto_pints(dim_hgto_bra,dim_hgto_ket, & 77 num_elec,num_mom,num_geo_pot) 78!f2py intent(in) :: orders_hgto_bra 79!f2py intent(in) :: coord_bra 80!f2py intent(in) :: exponent_bra 81!f2py intent(in) :: orders_hgto_ket 82!f2py intent(in) :: coord_ket 83!f2py intent(in) :: exponent_ket 84!f2py intent(in) :: order_elec 85!f2py intent(in) :: nucpot_origin 86!f2py intent(in) :: dipole_origin 87!f2py intent(in) :: scal_const 88!f2py intent(in) :: order_mom 89!f2py intent(in) :: order_geo_pot 90!f2py intent(in) :: dim_hgto_bra 91!f2py intent(in) :: dim_hgto_ket 92!f2py intent(in) :: num_elec 93!f2py intent(in) :: num_mom 94!f2py intent(in) :: num_geo_pot 95!f2py intent(out) :: hgto_pints 96!f2py depend(dim_hgto_bra) :: hgto_pints 97!f2py depend(dim_hgto_ket) :: hgto_pints 98!f2py depend(num_elec) :: hgto_pints 99!f2py depend(num_mom) :: hgto_pints 100!f2py depend(num_geo_pot) :: hgto_pints 101 integer orders_hket_elec(2) !orders of HGTOs on ket center with electronic derivatives 102 integer orders_hbra_mom(2) !orders of HGTOs on bra center with Cartesian multipole moments 103 integer max_order_nucpot !maximum order of geometric derivatives for \fn(nucpot_geom) 104 integer dim_nucpot !dimension of geometric derivatives of nuclear potential 105 !origin with zeroth order HGTOs 106 real(REALK), allocatable :: nucpot_pints(:) 107 !geometric derivatives of potential origin with zeroth order HGTOs 108 integer min_hgto_ket !minimum order of HGTOs on ket center after recovering HGTOs on ket center 109 integer max_geo_hket !maximum order of geometric derivatives after recovering HGTOs on ket center 110 integer dim_hket !dimension of HGTOs on ket center after recovering HGTOs on ket center 111 integer dim_geo_hket !dimension of geometric derivatives after recovering HGTOs on ket center 112 real(REALK), allocatable :: hket_pints(:,:) 113 !integrals after recovering HGTOs on ket center 114 integer dim_hbra_zero !dimension of HGTOs on bra center with zeroth order Cartesian multipole moment 115 integer dim_hket_zero !dimension of HGTOs on ket center with zeroth order Cartesian multipole moment 116 real(REALK), allocatable :: hbra_pints(:,:,:) 117 !integrals after recovering HGTOs on bra center 118 real(REALK), allocatable :: hmom_pints(:,:,:,:) 119 !integrals after recovering Cartesian multipole moments 120 integer ierr !error information 121#if defined(XTIME) 122 real(REALK) curr_time !current CPU time 123 ! sets current CPU time 124 call xtimer_set(curr_time) 125#endif 126 ! sets the orders of HGTOs on ket center with electronic derivatives 127 orders_hket_elec(1) = orders_hgto_ket(1)+order_elec 128 orders_hket_elec(2) = orders_hgto_ket(2)+order_elec 129 ! sets the orders of HGTOs on bra center with Cartesian multipole moments 130 orders_hbra_mom(1) = max(0,orders_hgto_bra(1)-order_mom) 131 orders_hbra_mom(2) = orders_hgto_bra(2)+order_mom 132 ! sets the maximum order and dimension of geometric derivatives for \fn(nucpot_geom) 133 max_order_nucpot = orders_hbra_mom(2)+orders_hket_elec(2)+order_geo_pot 134 if (order_geo_pot>0) then 135 dim_nucpot = ((max_order_nucpot+1)*(max_order_nucpot+2)*(max_order_nucpot+3) & 136 - order_geo_pot*(order_geo_pot+1)*(order_geo_pot+2))/6 137 else 138 dim_nucpot = (max_order_nucpot+1)*(max_order_nucpot+2)*(max_order_nucpot+3)/6 139 end if 140 ! allocates memory for the integrals from \fn(nucpot_geom) 141 allocate(nucpot_pints(dim_nucpot), stat=ierr) 142 if (ierr/=0) & 143 call error_stop("prim_hgto_nucpot", "failed to allocate nucpot_pints", & 144 dim_nucpot) 145 ! transfers Boys functions to geometric derivatives of nuclear potential 146 ! origin with zeroth order HGTOs 147 call nucpot_geom(coord_bra, exponent_bra, coord_ket, exponent_ket, & 148 nucpot_origin, scal_const, & 149 (/order_geo_pot,max_order_nucpot/), order_elec, & 150 dim_nucpot, nucpot_pints) 151 ! allocates memory for the integrals from \fn(nucpot_hket) 152 min_hgto_ket = max(orders_hket_elec(1)-orders_hbra_mom(2),0) 153 if (min_hgto_ket>0) then 154 dim_hket = ((orders_hket_elec(2)+1)*(orders_hket_elec(2)+2) & 155 *(orders_hket_elec(2)+3) & 156 - min_hgto_ket*(min_hgto_ket+1)*(min_hgto_ket+2))/6 157 else 158 dim_hket = (orders_hket_elec(2)+1)*(orders_hket_elec(2)+2) & 159 * (orders_hket_elec(2)+3)/6 160 end if 161 max_geo_hket = orders_hbra_mom(2)+order_geo_pot 162 if (order_geo_pot>0) then 163 dim_geo_hket = ((max_geo_hket+1)*(max_geo_hket+2)*(max_geo_hket+3) & 164 - order_geo_pot*(order_geo_pot+1)*(order_geo_pot+2))/6 165 else 166 dim_geo_hket = (max_geo_hket+1)*(max_geo_hket+2)*(max_geo_hket+3)/6 167 end if 168 allocate(hket_pints(dim_hket,dim_geo_hket), stat=ierr) 169 if (ierr/=0) & 170 call error_stop("prim_hgto_nucpot", "failed to allocate hket_pints", & 171 dim_hket*dim_geo_hket) 172 ! recovers the HGTOs on ket center 173 call nucpot_hket((/min_hgto_ket,orders_hket_elec(2)/), & 174 (/order_geo_pot,max_geo_hket/), & 175 coord_bra, exponent_bra, & 176 coord_ket, exponent_ket, & 177 dim_nucpot, nucpot_pints, & 178 dim_hket, dim_geo_hket, hket_pints) 179 deallocate(nucpot_pints) 180 if (order_mom==0) then 181 ! pure nuclear attraction 182 if (order_elec==0) then 183 ! recovers the HGTOs on bra center 184 call nucpot_hbra(orders_hbra_mom, orders_hket_elec, & 185 (/order_geo_pot,order_geo_pot/), & 186 coord_bra, exponent_bra, & 187 coord_ket, exponent_ket, & 188 dim_hket, dim_geo_hket, hket_pints, & 189 dim_hgto_bra, dim_hgto_ket, num_geo_pot, hgto_pints) 190 deallocate(hket_pints) 191 ! nuclear attraction + electronic derivatives 192 else 193 ! allocates the intermediate integrals from \fn(nucpot_hbra) 194 if (orders_hket_elec(1)==0) then 195 dim_hket_zero = (orders_hket_elec(2)+1)*(orders_hket_elec(2)+2) & 196 * (orders_hket_elec(2)+3)/6 197 else 198 dim_hket_zero = ((orders_hket_elec(2)+1)*(orders_hket_elec(2)+2) & 199 *(orders_hket_elec(2)+3) & 200 - orders_hket_elec(1)*(orders_hket_elec(1)+1) & 201 *(orders_hket_elec(1)+2))/6 202 end if 203 allocate(hbra_pints(dim_hgto_bra,dim_hket_zero,num_geo_pot), stat=ierr) 204 if (ierr/=0) & 205 call error_stop("prim_hgto_nucpot", "failed to allocate hbra_pints/ne", & 206 dim_hgto_bra*dim_hket_zero*num_geo_pot) 207 ! recovers the HGTOs on bra center 208 call nucpot_hbra(orders_hbra_mom, orders_hket_elec, & 209 (/order_geo_pot,order_geo_pot/), & 210 coord_bra, exponent_bra, & 211 coord_ket, exponent_ket, & 212 dim_hket, dim_geo_hket, hket_pints, & 213 dim_hgto_bra, dim_hket_zero, num_geo_pot, hbra_pints) 214 deallocate(hket_pints) 215 ! recovers the electronic derivatives 216 call scatter_multi_inner(dim_hgto_bra, dim_hket_zero, 1, num_geo_pot, & 217 hbra_pints, orders_hgto_ket, order_elec, & 218 dim_hgto_ket, num_elec, hgto_pints) 219 deallocate(hbra_pints) 220 end if 221 else 222 ! nuclear attraction + Cartesian multipole moments 223 if (order_elec==0) then 224 ! allocates the intermediate integrals from \fn(nucpot_hbra) 225 if (orders_hbra_mom(1)==0) then 226 dim_hbra_zero = (orders_hbra_mom(2)+1)*(orders_hbra_mom(2)+2) & 227 * (orders_hbra_mom(2)+3)/6 228 else 229 dim_hbra_zero = ((orders_hbra_mom(2)+1)*(orders_hbra_mom(2)+2) & 230 *(orders_hbra_mom(2)+3) & 231 - orders_hbra_mom(1)*(orders_hbra_mom(1)+1) & 232 *(orders_hbra_mom(1)+2))/6 233 end if 234 allocate(hbra_pints(dim_hbra_zero,dim_hgto_ket,num_geo_pot), stat=ierr) 235 if (ierr/=0) & 236 call error_stop("prim_hgto_nucpot", "failed to allocate hbra_pints/nm", & 237 dim_hbra_zero*dim_hgto_ket*num_geo_pot) 238 ! recovers the HGTOs on bra center 239 call nucpot_hbra(orders_hbra_mom, orders_hket_elec, & 240 (/order_geo_pot,order_geo_pot/), & 241 coord_bra, exponent_bra, & 242 coord_ket, exponent_ket, & 243 dim_hket, dim_geo_hket, hket_pints, & 244 dim_hbra_zero, dim_hgto_ket, num_geo_pot, hbra_pints) 245 deallocate(hket_pints) 246 ! recovers the Cartesian multipole moments 247 call london_mom_hgto(orders_hgto_bra, (/order_mom,order_mom/), & 248 coord_bra, exponent_bra, dipole_origin, & 249 1, dim_hbra_zero, dim_hgto_ket, num_geo_pot, & 250 hbra_pints, dim_hgto_bra, num_mom, hgto_pints) 251 deallocate(hbra_pints) 252 ! nuclear attraction + Cartesian multipole moments + electronic derivatives 253 else 254 ! allocates the intermediate integrals from \fn(nucpot_hbra) 255 if (orders_hbra_mom(1)==0) then 256 dim_hbra_zero = (orders_hbra_mom(2)+1)*(orders_hbra_mom(2)+2) & 257 * (orders_hbra_mom(2)+3)/6 258 else 259 dim_hbra_zero = ((orders_hbra_mom(2)+1)*(orders_hbra_mom(2)+2) & 260 *(orders_hbra_mom(2)+3) & 261 - orders_hbra_mom(1)*(orders_hbra_mom(1)+1) & 262 *(orders_hbra_mom(1)+2))/6 263 end if 264 if (orders_hket_elec(1)==0) then 265 dim_hket_zero = (orders_hket_elec(2)+1)*(orders_hket_elec(2)+2) & 266 * (orders_hket_elec(2)+3)/6 267 else 268 dim_hket_zero = ((orders_hket_elec(2)+1)*(orders_hket_elec(2)+2) & 269 *(orders_hket_elec(2)+3) & 270 - orders_hket_elec(1)*(orders_hket_elec(1)+1) & 271 *(orders_hket_elec(1)+2))/6 272 end if 273 allocate(hbra_pints(dim_hbra_zero,dim_hket_zero,num_geo_pot), stat=ierr) 274 if (ierr/=0) & 275 call error_stop("prim_hgto_nucpot", "failed to allocate hbra_pints/nme", & 276 dim_hbra_zero*dim_hket_zero*num_geo_pot) 277 ! recovers the HGTOs on bra center 278 call nucpot_hbra(orders_hbra_mom, orders_hket_elec, & 279 (/order_geo_pot,order_geo_pot/), & 280 coord_bra, exponent_bra, & 281 coord_ket, exponent_ket, & 282 dim_hket, dim_geo_hket, hket_pints, & 283 dim_hbra_zero, dim_hket_zero, num_geo_pot, hbra_pints) 284 deallocate(hket_pints) 285 ! recovers the Cartesian multipole moments 286 allocate(hmom_pints(dim_hgto_bra,dim_hket_zero,num_mom,num_geo_pot), stat=ierr) 287 if (ierr/=0) & 288 call error_stop("prim_hgto_nucpot", "failed to allocate hmom_pints", & 289 dim_hgto_bra*dim_hket_zero*num_mom*num_geo_pot) 290 call london_mom_hgto(orders_hgto_bra, (/order_mom,order_mom/), & 291 coord_bra, exponent_bra, dipole_origin, & 292 1, dim_hbra_zero, dim_hket_zero, num_geo_pot, & 293 hbra_pints, dim_hgto_bra, num_mom, hmom_pints) 294 deallocate(hbra_pints) 295 ! recovers the electronic derivatives 296 call scatter_multi_inner(dim_hgto_bra, dim_hket_zero, 1, num_mom*num_geo_pot, & 297 hmom_pints, orders_hgto_ket, order_elec , & 298 dim_hgto_ket, num_elec, hgto_pints) 299 deallocate(hmom_pints) 300 end if 301 end if 302#if defined(XTIME) 303 ! prints the CPU elapsed time 304 call xtimer_view(curr_time, "prim_hgto_nucpot", STDOUT) 305#endif 306 return 307 end subroutine prim_hgto_nucpot 308