1! This file is part of xtb.
2!
3! Copyright (C) 2019-2020 Sebastian Ehlert
4!
5! xtb is free software: you can redistribute it and/or modify it under
6! the terms of the GNU Lesser General Public License as published by
7! the Free Software Foundation, either version 3 of the License, or
8! (at your option) any later version.
9!
10! xtb is distributed in the hope that it will be useful,
11! but WITHOUT ANY WARRANTY; without even the implied warranty of
12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13! GNU Lesser General Public License for more details.
14!
15! You should have received a copy of the GNU Lesser General Public License
16! along with xtb.  If not, see <https://www.gnu.org/licenses/>.
17
18!> TODO
19module xtb_xtb_halogen
20   use xtb_mctc_accuracy, only : wp
21   use xtb_mctc_convert, only : aatoau
22   use xtb_lin, only : lin
23   use xtb_xtb_data
24   implicit none
25   private
26
27   public :: xbpot
28
29
30contains
31
32
33subroutine xbpot(halData,n,at,xyz,xblist,nxb,a,exb,g)
34   type(THalogenData), intent(in) :: halData
35   integer, intent(in) :: n
36   integer, intent(in) :: at(:)
37   integer, intent(in) :: nxb
38   integer, intent(in) :: xblist(:,:)
39   real(wp), intent(in) :: xyz(:,:)
40   real(wp), intent(inout) :: g(:,:)
41   real(wp), intent(inout) :: exb
42   real(wp), intent(in) :: a
43
44   integer :: m,k,AA,B,X,ati,atj
45   real(wp) :: cc,r0ax,t13,t14,t16
46   real(wp) :: d2ax,rax,term,aterm,xy,d2bx,d2ab,alp,lj2
47   real(wp) :: er,el,step,dxa(3),dxb(3),dba(3),dcosterm
48   real(wp) :: dtermlj,termlj,prefactor,numerator,denominator,rbx
49   alp=6.0_wp
50   lj2=0.50_wp*a
51
52   exb = 0.0_wp
53   if(nxb.lt.1) return
54
55   ! B-X...A
56   !$omp parallel do schedule(runtime) default(none) reduction(+:exb) &
57   !$omp shared(nxb, xblist, at, halData, xyz, lj2, alp, a) &
58   !$omp private(X, AA, B, ati, atj, cc, r0ax, dxa, dxb, dba, d2ax, &
59   !$omp& d2bx, d2ab, rax, XY, aterm, t13, t14, TERM)
60   do k=1,nxb
61      X =xblist(1,k)
62      AA=xblist(2,k)
63      B=xblist(3,k)
64      ati=at(X)
65      atj=at(AA)
66      cc=halData%bondStrength(ati)
67      r0ax=halData%radScale*(halData%atomicRad(ati)+halData%atomicRad(atj))
68      dxa=xyz(:,AA)-xyz(:,X)   ! acceptor - halogen
69      dxb=xyz(:, B)-xyz(:,X)   ! neighbor - halogen
70      dba=xyz(:,AA)-xyz(:, B)  ! acceptor - neighbor
71      d2ax=sum(dxa*dxa)
72      d2bx=sum(dxb*dxb)
73      d2ab=sum(dba*dba)
74      rax=sqrt(d2ax)
75      ! angle part. term = cos angle B-X-A
76      XY = SQRT(D2BX*D2AX)
77      TERM = (D2BX+D2AX-D2AB) / XY
78      aterm = (0.5_wp-0.25_wp*term)**alp
79      t13 = r0ax/rax
80      t14 = t13**a
81      exb = exb +  aterm*cc*(t14-halData%dampingPar*t13**lj2) / (1.0_wp+t14)
82   enddo
83
84   ! analytic gradient
85   !$omp parallel do schedule(runtime) default(none) reduction(+:g) &
86   !$omp shared(nxb, xblist, at, halData, xyz, lj2, alp) &
87   !$omp private(X, AA, B, ati, atj, cc, r0ax, dxa, dxb, dba, d2ax, &
88   !$omp& d2bx, d2ab, rax, XY, aterm, numerator, denominator, termLJ, &
89   !$omp& dtermlj, prefactor, dcosterm, t13, t14, rbx, TERM)
90   do k=1,nxb
91      X =xblist(1,k)
92      AA=xblist(2,k)
93      B=xblist(3,k)
94      ati=at(X)
95      atj=at(AA)
96      cc=halData%bondStrength(ati)
97      r0ax=halData%radScale*(halData%atomicRad(ati)+halData%atomicRad(atj))
98
99      dxa=xyz(:,AA)-xyz(:,X)   ! acceptor - halogen
100      dxb=xyz(:, B)-xyz(:,X)   ! neighbor - halogen
101      dba=xyz(:,AA)-xyz(:, B)  ! acceptor - neighbor
102
103      d2ax=sum(dxa*dxa)
104      d2bx=sum(dxb*dxb)
105      d2ab=sum(dba*dba)
106      rax=sqrt(d2ax)+1.0e-18_wp
107      rbx=sqrt(d2bx)+1.0e-18_wp
108
109      XY = SQRT(D2BX*D2AX)
110      TERM = (D2BX+D2AX-D2AB) / XY
111      ! now compute angular damping function
112      aterm = (0.5_wp-0.25_wp*term)**alp
113
114      ! set up weighted inverted distance and compute the modified Lennard-Jones potential
115      t14 = (r0ax/rax)**lj2 ! (rov/r)^lj2 ; lj2 = 6 in GFN1
116      numerator = (t14*t14 - halData%dampingPar*t14)
117      denominator = (1.0_wp + t14*t14)
118      termLJ= numerator/denominator
119
120      ! ----
121      ! LJ derivative
122      ! denominator part
123      dtermlj=2.0_wp*lj2*numerator*t14*t14/(rax*denominator*denominator)
124      ! numerator part
125      dtermlj=dtermlj+lj2*t14*(halData%dampingPar - 2.0_wp*t14)/(rax*denominator)
126      ! scale w/ angular damping term
127      dtermlj=dtermlj*aterm*cc/rax
128      ! gradient for the acceptor
129      g(:,AA)=g(:,AA)+dtermlj*dxa(:)
130      ! halogen gradient
131      g(:,X)=g(:,X)-dtermlj*dxa(:)
132      ! ----
133      ! cosine term derivative
134      prefactor=-0.250_wp*alp*(0.5_wp-0.25_wp*term)**(alp-1.0_wp)
135      prefactor=prefactor*cc*termlj
136      ! AX part
137      dcosterm=2.0_wp/rbx - term/rax
138      dcosterm=dcosterm*prefactor/rax
139      ! gradient for the acceptor
140      g(:,AA)=g(:,AA)+dcosterm*dxa(:)
141      ! halogen gradient
142      g(:,X)=g(:,X)-dcosterm*dxa(:)
143      ! BX part
144      dcosterm=2.0_wp/rax - term/rbx
145      dcosterm=dcosterm*prefactor/rbx
146      ! gradient for the acceptor
147      g(:,B)=g(:,B)+dcosterm*dxb(:)
148      ! halogen gradient
149      g(:,X)=g(:,X)-dcosterm*dxb(:)
150      ! AB part
151      t13=2.0_wp*prefactor/xy
152      ! acceptor
153      g(:,AA)=g(:,AA)-t13*dba(:)
154      ! neighbor
155      g(:,B)=g(:,B)+t13*dba(:)
156
157   enddo
158end subroutine xbpot
159
160
161end module xtb_xtb_halogen
162