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 tests spherical multipole integrals from Fortran 90 modules using SGTOs.
18!!
19!!  2012-03-20, Bin Gao
20!!  * first version
21
22  !> \brief tests spherical multipole integrals from Fortran 90 modules using SGTOs
23  !> \author Bin Gao
24  !> \date 2012-03-20
25  !> \param io_log is the IO unit of log file
26  !> \return test_failed indicates if the test is failed
27  subroutine test_f90mod_sgto_sphmom(io_log, test_failed)
28    use xkind
29    ! Fortran 90 module of Gen1Int library
30    use gen1int
31    ! module of HTML test log routines
32    use html_log
33    implicit none
34    integer, intent(in) :: io_log
35    logical, intent(out) :: test_failed
36! water molecule using cc-pV6Z basis sets
37#include "test_f90/water_ccpv6z.h"
38! test parameters and referenced results
39#include "water_f90mod_sgto_sphmom.h"
40    integer itst                        !incremental recorder over test cases
41    integer iblock, jblock              !incremental recorders over blocks of sub-shells
42    integer idx_bra, idx_ket            !indices of bra and ket centers
43    integer num_sgto_bra, num_sgto_ket  !number of spherical Gaussians of bra and ket centers
44    type(one_prop_t) one_prop           !one-electron property integral operator
45    integer num_prop                    !number of property integral matrices
46    real(REALK), allocatable :: contr_ints(:,:,:,:,:)
47                                        !contracted integrals
48    integer ierr                        !error information
49    integer addr_ref                    !address of referenced results
50    real(REALK) begin_time              !begin of CPU time
51    real(REALK) end_time                !end of CPU time
52    logical different                   !if result from \fn(contr_sgto_carmom) is different from reference
53    integer iopt                        !incremental recorder over operators
54    integer icontr, jcontr              !incremental recorders over contractions
55    integer isgto, jsgto                !incremental recorders over GTOs
56    ! assumes the test will pass
57    test_failed = .false.
58    call html_log_text("Tests spherical multipole integrals of water molecule "// &
59                       "using cc-pV6Z spherical GTOs ...", io_log)
60    ! initializes the information of spherical multipole integrals
61    call OnePropCreate(prop_name=INT_SPHER_MULTIPOLE, &
62                       one_prop=one_prop,             &
63                       info_prop=ierr,                &
64                       dipole_origin=DIP_ORIGIN,      &
65                       order_mom=ORDER_MOM)
66    ! gets the number of property integral matrices
67    call OnePropGetNumProp(one_prop=one_prop, num_prop=num_prop)
68    ! loops over test cases
69    addr_ref = 0
70    do itst = 1, NUM_TESTS
71      ! gets the blocks of sub-shells
72      iblock = BRA_BLOCK(itst)
73      jblock = KET_BLOCK(itst)
74      call html_log_int_number("ID of block of sub-shells on bra center:", &
75                               iblock, "I4", io_log)
76      call html_log_int_number("ID of block of sub-shells on ket center:", &
77                               jblock, "I4", io_log)
78      ! gets the begin time
79      call xtimer_set(begin_time)
80      ! gets the indices of bra and ket centers
81      idx_bra = IDX_CENT_BLOCK(iblock)
82      idx_ket = IDX_CENT_BLOCK(jblock)
83      ! sets the number of spherical Gaussians of bra and ket centers
84      num_sgto_bra = 2*ANGULAR_BLOCK(iblock)+1
85      num_sgto_ket = 2*ANGULAR_BLOCK(jblock)+1
86      ! allocates the contracted integrals
87      allocate(contr_ints(num_sgto_bra,NUM_CONTR_BLOCK(iblock), &
88                          num_sgto_ket,NUM_CONTR_BLOCK(jblock), &
89                          num_prop), stat=ierr)
90      if (ierr/=0) then
91        call html_log_int_number(log_text="Failed to allocate contr_ints:",    &
92                                 log_int=num_sgto_bra*NUM_CONTR_BLOCK(iblock)  &
93                                         *num_sgto_ket*NUM_CONTR_BLOCK(jblock) &
94                                         *num_prop,                            &
95                                 fmt_int="I12", io_log=io_log, font_color="red")
96        test_failed = .true.
97        cycle
98      end if
99      ! computes the contracted integrals
100      call OnePropGetIntegral(idx_bra=idx_bra, coord_bra=COORDS(3*idx_bra-2:3*idx_bra),  &
101                              angular_bra=ANGULAR_BLOCK(iblock),                         &
102                              num_prim_bra=NUM_PRIM_BLOCK(iblock),                       &
103                              exponent_bra=EXPONENT_BLOCK(START_PRIM_BLOCK(iblock)+1     &
104                                                          :START_PRIM_BLOCK(iblock+1)),  &
105                              num_contr_bra=NUM_CONTR_BLOCK(iblock),                     &
106                              contr_coef_bra=REF_NRM_SGTO(START_CONTR_BLOCK(iblock)+1    &
107                                                          :START_CONTR_BLOCK(iblock+1)), &
108                              idx_ket=idx_ket, coord_ket=COORDS(3*idx_ket-2:3*idx_ket),  &
109                              angular_ket=ANGULAR_BLOCK(jblock),                         &
110                              num_prim_ket=NUM_PRIM_BLOCK(jblock),                       &
111                              exponent_ket=EXPONENT_BLOCK(START_PRIM_BLOCK(jblock)+1     &
112                                                          :START_PRIM_BLOCK(jblock+1)),  &
113                              num_contr_ket=NUM_CONTR_BLOCK(jblock),                     &
114                              contr_coef_ket=REF_NRM_SGTO(START_CONTR_BLOCK(jblock)+1    &
115                                                          :START_CONTR_BLOCK(jblock+1)), &
116                              spher_gto=.true., one_prop=one_prop,                       &
117                              num_gto_bra=num_sgto_bra, num_gto_ket=num_sgto_ket,        &
118                              num_opt=num_prop, contr_ints=contr_ints)
119      ! gets the end time
120      call xtimer_set(end_time)
121      ! prints the CPU elapsed time
122      call html_log_real_number(log_text="Time (s) used for OnePropGetIntegral:", &
123                                log_real=end_time-begin_time, fmt_real="F10.4",   &
124                                io_log=io_log, font_color="blue")
125      ! checks the results
126      do iopt = 1, num_prop
127        do jcontr = 1, NUM_CONTR_BLOCK(jblock)
128          do jsgto = 1, num_sgto_ket
129            do icontr = 1, NUM_CONTR_BLOCK(iblock)
130              do isgto = 1, num_sgto_bra
131                addr_ref = addr_ref+1
132                call check_difference(REF_CONTR_INTS(addr_ref),                   &
133                                      contr_ints(isgto,icontr,jsgto,jcontr,iopt), &
134                                      different)
135                if (different) then
136                  call html_log_int_number("ID of test:", itst, "I3", io_log)
137                  call html_log_int_array(log_text="ID of contraction and operator:", &
138                                          log_int=(/isgto,icontr,jsgto,jcontr,iopt/), &
139                                          fmt_int="I3", io_log=io_log)
140                  call html_log_real_number(                                  &
141                         log_text="Referenced spherical multipole integral:", &
142                         log_real=REF_CONTR_INTS(addr_ref),                   &
143                         fmt_real="Es16.8", io_log=io_log, font_color="blue")
144                  call html_log_real_number(                                  &
145                         log_text="Result from OnePropGetIntegral:",          &
146                         log_real=contr_ints(isgto,icontr,jsgto,jcontr,iopt), &
147                         fmt_real="Es16.8", io_log=io_log, font_color="red")
148                  test_failed = .true.
149                end if
150              end do
151            end do
152          end do
153        end do
154      end do
155      ! cleans
156      deallocate(contr_ints)
157    end do
158    return
159  end subroutine test_f90mod_sgto_sphmom
160