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!!  Tests subroutines in contr_cgto_carmom.F90
18!!
19!!  2011-04-03, Bin Gao
20!!  * first version
21
22  !> \brief tests subroutines in contr_cgto_carmom.F90
23  !> \author Bin Gao
24  !> \date 2011-04-03
25  !> \param io_log is the IO unit of log file
26  !> \return test_failed indicates if the test is failed
27  subroutine test_contr_cgto_carmom(io_log, test_failed)
28    use xkind
29    ! module of HTML test log routines
30    use html_log
31    implicit none
32    integer, intent(in) :: io_log
33    logical, intent(out) :: test_failed
34! water molecule using cc-pV6Z basis sets
35#include "test_f90/water_ccpv6z.h"
36! test parameters and referenced results
37#include "water_ccpv6z_cgto_carmom.h"
38    integer itst, jtst, ktst, ltst      !incremental recorders over test cases
39    integer iblock, jblock              !incremental recorders over blocks of sub-shells
40    integer num_opt                     !number of operators including derivatives
41    real(REALK), allocatable :: contr_ints(:,:,:,:,:)
42                                        !contracted integrals
43    integer idx_bra, idx_ket            !indices of bra and ket centers
44    integer num_cgto_bra, num_cgto_ket  !number of Cartesian Gaussians of bra and ket centers
45    integer ierr                        !error information
46    integer addr_ref                    !address of referenced results
47    real(REALK) begin_time              !begin of CPU time
48    real(REALK) end_time                !end of CPU time
49    logical different                   !if result from \fn(contr_cgto_carmom) is different from reference
50    real(REALK), allocatable :: ro_cints(:,:,:,:,:)  !reordered integrals
51    integer, allocatable :: power_bra(:,:)           !powers on bra and ket centers
52    integer, allocatable :: power_ket(:,:)
53    integer xpow, ypow                               !incremental recorders over xy powers
54    integer addr_bra, addr_ket                       !addresses of GTOs on bra and ket centers
55    integer icontr, jcontr, icgto, jcgto, iopt       !incremental recorders
56    ! we first test subroutines in reorder_ints.F90 which will be used later on
57    call html_log_text("We first test subroutines in reorder_ints.F90 ...", io_log)
58    call test_reorder_ints(io_log, test_failed)
59    if (test_failed) then
60      call html_log_text(log_text="Tests subroutines in reorder_ints.F90 failed! "// &
61                                  "No further test will be performed!",              &
62                         io_log=io_log, font_color="red")
63      return
64    else
65      call html_log_text("Tests subroutines in reorder_ints.F90 passed!", io_log)
66    end if
67    ! dumps information of tests
68    call html_log_text("Tests water molecule using cc-pV6Z Cartesian GTOs ...", io_log)
69    ! loops over test cases
70    jtst = 1
71    addr_ref = 0
72    do itst = 1, NUM_TESTS
73      ktst = jtst+2
74      call html_log_int_number("Order of electronic derivatives:", &
75                               ORDER_ELEC(itst), "I4", io_log)
76      if (IDX_DIPORG(itst)>0)                               &
77        call html_log_int_number("Index of dipole origin:", &
78                                 IDX_DIPORG(itst), "I4", io_log)
79      call html_log_real_array("Coordinates of dipole origin:", DIP_ORIGIN(jtst:ktst), &
80                               "F16.8", io_log)
81      call html_log_real_number("Scaling constant:", SCAL_CONST(itst), "F16.8", io_log)
82      call html_log_int_number("Order of Cartesian multipole moments:", &
83                               ORDER_MOM(itst), "I4", io_log)
84      call html_log_int_number("Order of partial geometric derivatives on bra center:", &
85                               ORDER_GEO_BRA(itst), "I4", io_log)
86      call html_log_int_number("Order of partial geometric derivatives on ket center:", &
87                               ORDER_GEO_KET(itst), "I4", io_log)
88      call html_log_int_number("Order of geometric derivatives on dipole origin:", &
89                               ORDER_GEO_MOM(itst), "I4", io_log)
90      ! computes the number of total geometric derivatives
91      num_opt = 1
92      do ltst = jtst, jtst+NUM_CENTS(itst)-1
93        num_opt = num_opt*(ORDER_CENT(ltst)+1)*(ORDER_CENT(ltst)+2)/2
94      end do
95      ltst = jtst+NUM_CENTS(itst)-1
96      if (NUM_CENTS(itst)>0) then
97        call html_log_int_number("Numbers of differentiated centers of total "// &
98                                 "geometric derivatives:", NUM_CENTS(itst), "I4", io_log)
99        call html_log_int_array("Indices of differentiated centers:", &
100                                IDX_CENT(jtst:ltst), "I4", io_log)
101        call html_log_int_array("Orders of differentiated centers:", &
102                                ORDER_CENT(jtst:ltst), "I4", io_log)
103        if (NUM_CENTS(itst)>3) then
104          call html_log_int_number(log_text="Invalid number of differentiated "//      &
105                                            "centers of total geometric derivatives:", &
106                                   log_int=NUM_CENTS(itst), fmt_int="I4",              &
107                                   io_log=io_log, font_color="red")
108          test_failed = .true.
109          cycle
110        end if
111      end if
112      ! gets the begin time
113      call xtimer_set(begin_time)
114      ! computs the number of operators including the derivatives
115      num_opt = num_opt                                         &
116              * (ORDER_ELEC(itst)+1)*(ORDER_ELEC(itst)+2)       &
117              * (ORDER_MOM(itst)+1)*(ORDER_MOM(itst)+2)         &
118              * (ORDER_GEO_BRA(itst)+1)*(ORDER_GEO_BRA(itst)+2) &
119              * (ORDER_GEO_KET(itst)+1)*(ORDER_GEO_KET(itst)+2) &
120              * (ORDER_GEO_MOM(itst)+1)*(ORDER_GEO_MOM(itst)+2)/32
121      ! gets the blocks of sub-shells
122      iblock = BRA_BLOCK(itst)
123      jblock = KET_BLOCK(itst)
124      call html_log_int_number("ID of block of sub-shells on bra center:", &
125                               iblock, "I4", io_log)
126      call html_log_int_number("ID of block of sub-shells on ket center:", &
127                               jblock, "I4", io_log)
128      ! gets the indices of bra and ket centers
129      idx_bra = IDX_CENT_BLOCK(iblock)
130      idx_ket = IDX_CENT_BLOCK(jblock)
131      ! sets the number of Cartesian Gaussians of bra and ket centers
132      num_cgto_bra = (ANGULAR_BLOCK(iblock)+1)*(ANGULAR_BLOCK(iblock)+2)/2
133      num_cgto_ket = (ANGULAR_BLOCK(jblock)+1)*(ANGULAR_BLOCK(jblock)+2)/2
134      ! allocates the contracted integrals
135      allocate(contr_ints(num_cgto_bra,NUM_CONTR_BLOCK(iblock), &
136                          num_cgto_ket,NUM_CONTR_BLOCK(jblock), &
137                          num_opt), stat=ierr)
138      if (ierr/=0) then
139        call html_log_int_number(log_text="Failed to allocate contr_ints:",    &
140                                 log_int=num_cgto_bra*NUM_CONTR_BLOCK(iblock)  &
141                                         *num_cgto_ket*NUM_CONTR_BLOCK(jblock) &
142                                         *num_opt,                             &
143                                 fmt_int="I12", io_log=io_log, font_color="red")
144        test_failed = .true.
145        cycle
146      end if
147      ! computes the contracted integrals
148      call contr_cgto_carmom(idx_bra, COORDS(3*idx_bra-2:3*idx_bra),        &
149                             ANGULAR_BLOCK(iblock), NUM_PRIM_BLOCK(iblock), &
150                             EXPONENT_BLOCK(START_PRIM_BLOCK(iblock)+1      &
151                                            :START_PRIM_BLOCK(iblock+1)),   &
152                             NUM_CONTR_BLOCK(iblock),                       &
153                             REF_NRM_CGTO(START_CONTR_BLOCK(iblock)+1       &
154                                          :START_CONTR_BLOCK(iblock+1)),    &
155                             idx_ket, COORDS(3*idx_ket-2:3*idx_ket),        &
156                             ANGULAR_BLOCK(jblock), NUM_PRIM_BLOCK(jblock), &
157                             EXPONENT_BLOCK(START_PRIM_BLOCK(jblock)+1      &
158                                            :START_PRIM_BLOCK(jblock+1)),   &
159                             NUM_CONTR_BLOCK(jblock),                       &
160                             REF_NRM_CGTO(START_CONTR_BLOCK(jblock)+1       &
161                                          :START_CONTR_BLOCK(jblock+1)),    &
162                             ORDER_ELEC(itst), IDX_DIPORG(itst),            &
163                             DIP_ORIGIN(jtst:ktst), SCAL_CONST(itst),       &
164                             ORDER_MOM(itst), ORDER_GEO_BRA(itst),          &
165                             ORDER_GEO_KET(itst), ORDER_GEO_MOM(itst),      &
166                             NUM_CENTS(itst), IDX_CENT(jtst:ltst),          &
167                             ORDER_CENT(jtst:ltst), num_cgto_bra,           &
168                             num_cgto_ket, num_opt, contr_ints)
169      ! gets the end time
170      call xtimer_set(end_time)
171      ! prints the CPU elapsed time
172      call html_log_real_number(log_text="Time (s) used for contr_cgto_carmom:", &
173                                log_real=end_time-begin_time, fmt_real="F10.4",  &
174                                io_log=io_log, font_color="blue")
175      ! prepares the powers on bra and ket centers
176      allocate(power_bra(3,num_cgto_bra), stat=ierr)
177      if (ierr/=0) then
178        call html_log_int_number(log_text="Failed to allocate power_bra:", &
179                                 log_int=3*num_cgto_bra, fmt_int="I12",    &
180                                 io_log=io_log, font_color="red")
181        deallocate(contr_ints)
182        test_failed = .true.
183        cycle
184      end if
185      allocate(power_ket(3,num_cgto_ket), stat=ierr)
186      if (ierr/=0) then
187        call html_log_int_number(log_text="Failed to allocate power_ket:", &
188                                 log_int=3*num_cgto_ket, fmt_int="I12",    &
189                                 io_log=io_log, font_color="red")
190        deallocate(power_bra)
191        deallocate(contr_ints)
192        test_failed = .true.
193        cycle
194      end if
195      addr_bra = 0
196      do xpow = ANGULAR_BLOCK(iblock), 0, -1
197        do ypow = ANGULAR_BLOCK(iblock)-xpow, 0, -1
198          addr_bra = addr_bra+1
199          power_bra(1,addr_bra) = xpow
200          power_bra(2,addr_bra) = ypow
201          power_bra(3,addr_bra) = ANGULAR_BLOCK(iblock)-(xpow+ypow)
202        end do
203      end do
204      addr_ket = 0
205      do xpow = ANGULAR_BLOCK(jblock), 0, -1
206        do ypow = ANGULAR_BLOCK(jblock)-xpow, 0, -1
207          addr_ket = addr_ket+1
208          power_ket(1,addr_ket) = xpow
209          power_ket(2,addr_ket) = ypow
210          power_ket(3,addr_ket) = ANGULAR_BLOCK(jblock)-(xpow+ypow)
211        end do
212      end do
213      ! reorders the integrals in the order of Dalton format
214      allocate(ro_cints(num_cgto_bra,NUM_CONTR_BLOCK(iblock), &
215                        num_cgto_ket,NUM_CONTR_BLOCK(jblock), &
216                        num_opt), stat=ierr)
217      if (ierr/=0) then
218        call html_log_int_number(log_text="Failed to allocate ro_cints:",      &
219                                 log_int=num_cgto_bra*NUM_CONTR_BLOCK(iblock)  &
220                                         *num_cgto_ket*NUM_CONTR_BLOCK(jblock) &
221                                         *num_opt,                             &
222                                 fmt_int="I12", io_log=io_log, font_color="red")
223        deallocate(power_bra)
224        deallocate(power_ket)
225        deallocate(contr_ints)
226        test_failed = .true.
227        cycle
228      end if
229      call reorder_cgto_ints(ANGULAR_BLOCK(iblock), num_cgto_bra, power_bra,   &
230                             ANGULAR_BLOCK(jblock), num_cgto_ket, power_ket,   &
231                             NUM_CONTR_BLOCK(iblock), NUM_CONTR_BLOCK(jblock), &
232                             num_opt, contr_ints, ro_cints)
233      ! cleans
234      deallocate(contr_ints)
235      deallocate(power_bra)
236      deallocate(power_ket)
237      ! checks the results
238      do iopt = 1, num_opt
239        do jcontr = 1, NUM_CONTR_BLOCK(jblock)
240          do jcgto = 1, num_cgto_ket
241            do icontr = 1, NUM_CONTR_BLOCK(iblock)
242              do icgto = 1, num_cgto_bra
243                addr_ref = addr_ref+1
244                call check_difference(REF_CONTR_INTS(addr_ref),                 &
245                                      ro_cints(icgto,icontr,jcgto,jcontr,iopt), &
246                                      different)
247                if (different) then
248                  call html_log_int_number("ID of test:", itst, "I3", io_log)
249                  call html_log_int_array(log_text="ID of contraction and operator:", &
250                                          log_int=(/icgto,icontr,jcgto,jcontr,iopt/), &
251                                          fmt_int="I3", io_log=io_log)
252                  call html_log_real_number(                         &
253                         log_text="Referenced contracted integral:", &
254                         log_real=REF_CONTR_INTS(addr_ref),          &
255                         fmt_real="Es16.8", io_log=io_log, font_color="blue")
256                  call html_log_real_number(                                &
257                         log_text="Result from contr_cgto_carmom:",         &
258                         log_real=ro_cints(icgto,icontr,jcgto,jcontr,iopt), &
259                         fmt_real="Es16.8", io_log=io_log, font_color="red")
260                  test_failed = .true.
261                end if
262              end do
263            end do
264          end do
265        end do
266      end do
267      ! cleans
268      deallocate(ro_cints)
269      ! updates the incremental recorder
270      jtst = jtst+3
271    end do
272    return
273  end subroutine test_contr_cgto_carmom
274