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