1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Utilities for evaluating the residual part (1/r^3) of Integrals for
8!>        semi-empiric methods
9!> \author Teodoro Laino (11.2008) [tlaino]
10! **************************************************************************************************
11MODULE semi_empirical_int3_utils
12
13   USE input_constants,                 ONLY: do_method_pchg
14   USE kinds,                           ONLY: dp
15   USE semi_empirical_int_arrays,       ONLY: clm_d,&
16                                              indexb
17   USE semi_empirical_types,            ONLY: semi_empirical_type
18#include "./base/base_uses.f90"
19
20   IMPLICIT NONE
21   PRIVATE
22   LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
23   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int3_utils'
24
25   PUBLIC ::   ijkl_low_3, charg_int_3, dcharg_int_3, coeff_int_3
26
27CONTAINS
28
29! **************************************************************************************************
30!> \brief Low level general driver for computing residual part of semi-empirical
31!>        integrals <ij|kl> and their derivatives
32!>        The residual part is the leading 1/r^3 term
33!>
34!> \param sepi ...
35!> \param sepj ...
36!> \param ij ...
37!> \param kl ...
38!> \param li ...
39!> \param lj ...
40!> \param lk ...
41!> \param ll ...
42!> \param ic ...
43!> \param itype ...
44!> \param eval ...
45!> \return ...
46!> \date 11.2008 [tlaino]
47!> \author Teodoro Laino [tlaino]
48! **************************************************************************************************
49   FUNCTION ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, itype, eval) RESULT(res)
50      TYPE(semi_empirical_type), POINTER                 :: sepi, sepj
51      INTEGER, INTENT(IN)                                :: ij, kl, li, lj, lk, ll, ic, itype
52      REAL(KIND=dp)                                      :: eval, res
53
54      CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_low_3', routineP = moduleN//':'//routineN
55
56      INTEGER                                            :: l1, l2, lij, lkl
57      REAL(KIND=dp)                                      :: add, ccc, chrg, pij, pkl, sum
58
59      sum = 0.0_dp
60      l1 = ABS(li - lj)
61      lij = indexb(li + 1, lj + 1)
62      l2 = ABS(lk - ll)
63      lkl = indexb(lk + 1, ll + 1)
64
65      ! Standard value of the integral
66      IF (l1 == 0) THEN
67         IF (lij == 1) THEN
68            pij = sepi%ko(1)
69            IF (ic == 1) THEN
70               pij = sepi%ko(9)
71            END IF
72         ELSE IF (lij == 3) THEN
73            pij = sepi%ko(7)
74         ELSE IF (lij == 6) THEN
75            pij = sepi%ko(8)
76         END IF
77      END IF
78      !
79      IF (l2 == 0) THEN
80         IF (lkl == 1) THEN
81            pkl = sepj%ko(1)
82            IF (ic == 2) THEN
83               pkl = sepj%ko(9)
84            END IF
85         ELSE IF (lkl == 3) THEN
86            pkl = sepj%ko(7)
87         ELSE IF (lkl == 6) THEN
88            pkl = sepj%ko(8)
89         END IF
90      END IF
91      IF (l1 == 0 .AND. l2 == 0) THEN
92         IF (itype == do_method_pchg) THEN
93            add = 0.0_dp
94         ELSE
95            add = (pij + pkl)**2
96         END IF
97         ccc = clm_d(ij, l1, 0)*clm_d(kl, l2, 0)
98         IF (ABS(ccc) > EPSILON(0.0_dp)) THEN
99            chrg = eval(l1, l2, add)
100            sum = chrg
101         END IF
102      END IF
103      res = sum
104   END FUNCTION ijkl_low_3
105
106! **************************************************************************************************
107!> \brief Evaluates the residual Interaction function between two point-charges
108!>        The term evaluated is the 1/r^3 (for short range interactions)
109!>        r    -  Distance r12
110!>        l1   -  Quantum numbers for multipole of configuration 1
111!>        l2   -  Quantum numbers for multipole of configuration 2
112!>        add  -  additive term
113!>
114!> \param r ...
115!> \param l1 ...
116!> \param l2 ...
117!> \param add ...
118!> \return ...
119!> \date 11.2008 [tlaino]
120!> \author Teodoro Laino [tlaino]
121! **************************************************************************************************
122   FUNCTION charg_int_3(r, l1, l2, add) RESULT(charg)
123      REAL(KIND=dp), INTENT(in)                          :: r
124      INTEGER, INTENT(in)                                :: l1, l2
125      REAL(KIND=dp), INTENT(in)                          :: add
126      REAL(KIND=dp)                                      :: charg
127
128      CHARACTER(len=*), PARAMETER :: routineN = 'charg_int_3', routineP = moduleN//':'//routineN
129
130! Computing only residual Integral Values
131
132      charg = 0.0_dp
133      ! Q - Q.
134      IF (l1 == 0 .AND. l2 == 0) THEN
135         charg = -add/(2.0_dp*r**3)
136         RETURN
137      END IF
138      ! We should NEVER reach this point
139      CPABORT("")
140   END FUNCTION charg_int_3
141
142! **************************************************************************************************
143!> \brief Evaluates the coefficient for the residual Interaction function
144!>        between two point-charges
145!>        l1   -  Quantum numbers for multipole of configuration 1
146!>        l2   -  Quantum numbers for multipole of configuration 2
147!>        add  -  additive term
148!>
149!> \param l1 ...
150!> \param l2 ...
151!> \param add ...
152!> \return ...
153!> \date 11.2008 [tlaino]
154!> \author Teodoro Laino [tlaino]
155! **************************************************************************************************
156   FUNCTION coeff_int_3(l1, l2, add) RESULT(coeff)
157      INTEGER, INTENT(in)                                :: l1, l2
158      REAL(KIND=dp), INTENT(in)                          :: add
159      REAL(KIND=dp)                                      :: coeff
160
161      CHARACTER(len=*), PARAMETER :: routineN = 'coeff_int_3', routineP = moduleN//':'//routineN
162
163! Computing only residual Integral Values
164
165      coeff = 0.0_dp
166      ! Q - Q.
167      IF (l1 == 0 .AND. l2 == 0) THEN
168         coeff = -add/2.0_dp
169         RETURN
170      END IF
171      ! We should NEVER reach this point
172      CPABORT("")
173   END FUNCTION coeff_int_3
174
175! **************************************************************************************************
176!> \brief Derivatives of residual interaction function between two point-charges
177!>
178!>        r    -  Distance r12
179!>        l1   -  Quantum numbers for multipole of configuration 1
180!>        l2   -  Quantum numbers for multipole of configuration 2
181!>        add  -  additive term
182!>
183!> \param r ...
184!> \param l1 ...
185!> \param l2 ...
186!> \param add ...
187!> \return ...
188!> \date 11.2008 [tlaino]
189!> \author Teodoro Laino [tlaino]
190! **************************************************************************************************
191   FUNCTION dcharg_int_3(r, l1, l2, add) RESULT(charg)
192      REAL(KIND=dp), INTENT(in)                          :: r
193      INTEGER, INTENT(in)                                :: l1, l2
194      REAL(KIND=dp), INTENT(in)                          :: add
195      REAL(KIND=dp)                                      :: charg
196
197      CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_3', routineP = moduleN//':'//routineN
198
199! Computing only residual Integral Derivatives
200
201      charg = 0.0_dp
202      ! Q - Q.
203      IF (l1 == 0 .AND. l2 == 0) THEN
204         charg = 3.0_dp*add/(2.0_dp*r**4)
205         RETURN
206      END IF
207      ! We should NEVER reach this point
208      CPABORT("")
209   END FUNCTION dcharg_int_3
210
211END MODULE semi_empirical_int3_utils
212