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