1! This file is part of xtb. 2! 3! Copyright (C) 2017-2020 Stefan Grimme 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 18module xtb_aespot 19 use xtb_mctc_accuracy, only : wp 20 use xtb_intpack, only : olap,divpt,rhftce,prod,opab1,opab4,propa 21 use xtb_xtb_data 22 integer,private, parameter :: llao (0:3) = (/ 1, 3, 6,10/) 23 integer,private, parameter :: llao2(0:3) = (/ 1, 3, 5, 7/) 24 25contains 26 27! setdqlist: precomputes the dipole and quadrupole potential terms for Fock matrix 28! thr : neglect matrix elements in potential 29! nao : # of spherical AOs (SAOs) 30! dpint : dipole integral matrix, dimension 3,nao*(nao+1)/2 31! qpint : quadrupole integral matrix, dimension 6,nao*(nao+1)/2 32! ndp,nqp : number of elements to be computed in Fock matrix with X-dip and X-qpole terms 33! matdlst,matqlst : index list, to which AO, the ndp/nqp potential terms refer to 34subroutine setdqlist(nao,ndp,nqp,thr,dpint,qpint,matdlst,matqlst) 35 implicit none 36 integer, intent(in) :: nao 37 integer, intent(inout) :: ndp,nqp 38 real(wp),intent(in) :: thr 39 real(wp),intent(in) :: dpint(3,nao,nao) 40 real(wp),intent(in) :: qpint(6,nao,nao) 41 integer, intent(inout):: matqlst(2,nqp),matdlst(2,ndp) 42 43 real(wp) tmp1,tmp2,tmp3,tmp4 44 real(wp) skj,r1,r2,tt,t1,t2,t3,t4,thr2,f 45 ! stuff for potential 46 47 integer i,j,k,l,m,ii,jj,ll,kk,mq,md,ij 48 49 ! INFO: this threshold must be slightly larger than max(0,thr2), 50 ! where thr2 is the one used in screening in routine aesdqint 51 thr2 = thr*1.0d-2 ! we compare squared int-elements 52 ! thr2 = 1.0d-20 ! conservative, keep all terms 53 54 md = 0 55 mq = 0 56 ! set uo matrix lists 57 ij = 0 58 do i = 1,nao 59 do j = 1,i 60 ij = ij+1 61 tmp1 = 0.0_wp 62 tmp2 = 0.0_wp 63 kk = 0 64 do k = 1,3 65 tmp1 = tmp1+dpint(k,i,j)*dpint(k,i,j) 66 tmp2 = tmp2-qpint(k,i,j)*qpint(k,i,j) 67 enddo 68 do k = 1,6 69 tmp2 = tmp2+2.0_wp*qpint(k,i,j)*qpint(k,i,j) 70 enddo 71 if(tmp1.gt.thr2)then 72 md = md+1 73 matdlst(1,md) = int(i,2) 74 matdlst(2,md) = int(j,2) 75 endif 76 if(tmp2.gt.thr2)then 77 mq = mq+1 78 matqlst(1,mq) = int(i,2) 79 matqlst(2,mq) = int(j,2) 80 endif 81 enddo 82 enddo 83 ndp = md 84 nqp = mq 85end subroutine setdqlist 86 87! scalecamm: scale all anisotropic CAMMs by element-specific parameters 88! nat : # of atoms 89! at(nat) : atom-to-element identifier 90! dipm(3,nat) : cumulative atomic dipole moments (x,y,z) 91! qp(6,nat) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 92subroutine scalecamm(nat,at,dipm,qp) 93 implicit none 94 integer, intent(in) :: nat,at(nat) 95 real(wp), intent(inout):: dipm(3,nat),qp(6,nat) 96 integer i,iat 97 98 ! CAMM scaling 99 do i = 1,nat 100 qp(1:6,i) = qp(1:6,i)*(3./3.) 101 enddo 102end subroutine scalecamm 103 104! unscalecamm: unscale all anisotropic CAMMs from element-specific parameters 105! nat : # of atoms 106! at(nat) : atom-to-element identifier 107! dipm(3,nat) : cumulative atomic dipole moments (x,y,z) 108! qp(6,nat) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 109subroutine unscalecamm(nat,at,dipm,qp) 110 implicit none 111 integer, intent(in) :: nat,at(nat) 112 real(wp), intent(inout):: dipm(3,nat),qp(6,nat) 113 real(wp) aesi 114 integer i,iat 115 116 ! CAMM scaling 117 do i = 1,nat 118 qp(1:6,i) = qp(1:6,i)*(3./3.) 119 enddo 120end subroutine unscalecamm 121 122 123 124 125! mmpop: compute the cumulative atomic dipole and quadrupole moments via Mulliken population analysis 126! nat : # of atoms 127! nao : # of spherical AOs (SAOs) 128! aoat2(nao) : SAO to atom intex 129! s(nao,nao) : overlap matrix 130! xyz(3,nat) : cartesian coordinates 131! dpint : dipole integral matrix, dimension 3,nao*(nao+1)/2 132! qpint : quadrupole integral matrix, dimension 6,nao*(nao+1)/2 133! p(nao,nao) : density matrix 134! dipm(3,nat) : cumulative atomic dipole moments (x,y,z) 135! qp(6,nat) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 136subroutine mmompop(nat,nao,aoat2,xyz,p,s,dpint,qpint,dipm,qp) 137 implicit none 138 integer, intent(in) :: nao,nat,aoat2(:) 139 real(wp), intent(in) :: s(:, :) 140 real(wp), intent(in) :: p(:, :) 141 real(wp), intent(in) :: dpint(:, :, :) 142 real(wp), intent(in) :: qpint(:, :, :) 143 real(wp), intent(in) :: xyz(:, :) 144 real(wp), intent(out):: dipm(:, :) 145 real(wp), intent(out):: qp(:, :) 146 147 real(wp) xk1,xl1,xk2,xl2,pij,tii,tjj 148 real(wp) pqm,pdmk,pdml,ps,ra(3) 149 150 integer i,j,k,l,ii,jj,kl,kj,lin 151 152 !$acc enter data create(dipm(:, :), qp(:, :)) 153 154 ! CAMM 155 !$acc kernels default(present) 156 dipm = 0.0_wp 157 qp = 0.0_wp 158 !$acc end kernels 159 160 !$acc enter data copyin(nao, nat, aoat2(:), s(:, :), p(:, :), dpint(:, :, :), & 161 !$acc& qpint(:, :, :),xyz(:, :)) 162 163 !$acc parallel private(pij,ps,ra,k,l,ii,jj) 164 165 !$acc loop gang vector collapse(2) 166 do i = 1,nao 167 do j = 1,nao 168 if (j >= i) cycle 169 ii = aoat2(i) 170 jj = aoat2(j) 171 ra(1:3) = xyz(1:3,ii) 172 pij = p(j,i) 173 ps = pij*s(j,i) 174 ! the qpint is stored as xx,yy,zz,xy,xz,yz (from integral routine) 175 ! when doing the Mulliken population, we switch to lin-compatible sorting 176 ! i,e. xx,xy,yy,xz,yz,zz 177 !$acc loop vector private(xk1,xl1,xk2,xl2,tii,tjj,pqm,pdmk,pdml,kl,kj) 178 do k = 1,3 179 xk1 = ra(k) 180 xk2 = xyz(k,jj) 181 pdmk = pij*dpint(k,j,i) 182 tii = xk1*ps-pdmk 183 tjj = xk2*ps-pdmk 184 !$acc atomic 185 dipm(k,jj) = dipm(k,jj)+tjj 186 !$acc atomic 187 dipm(k,ii) = dipm(k,ii)+tii 188 ! off-diagonal 189 do l = 1,k-1 190 kl = k*(k-1)/2+l 191 kj = k+l+1 192 xl1 = ra(l) 193 xl2 = xyz(l,jj) 194 pdml = pij*dpint(l,j,i) 195 pqm = pij*qpint(kj,j,i) 196 tii = pdmk*xl1+pdml*xk1-xl1*xk1*ps-pqm 197 tjj = pdmk*xl2+pdml*xk2-xl2*xk2*ps-pqm 198 !$acc atomic 199 qp(kl,jj) = qp(kl,jj)+tjj 200 !$acc atomic 201 qp(kl,ii) = qp(kl,ii)+tii 202 enddo 203 ! diagonal 204 kl = k*(k+1)/2 205 pqm = pij*qpint(k,j,i) 206 tii = 2.0_wp*pdmk*xk1-xk1*xk1*ps-pqm 207 tjj = 2.0_wp*pdmk*xk2-xk2*xk2*ps-pqm 208 !$acc atomic 209 qp(kl,jj) = qp(kl,jj)+tjj 210 !$acc atomic 211 qp(kl,ii) = qp(kl,ii)+tii 212 enddo 213 enddo 214 enddo 215 216 !$acc loop gang vector 217 do i = 1,nao 218 ii = aoat2(i) 219 ra(1:3) = xyz(1:3,ii) 220 pij = p(i,i) 221 ps = pij*s(i,i) 222 ! the qpint is stored as xx,yy,zz,xy,xz,yz (from integral routine) 223 ! when doing the Mulliken population, we switch to lin-compatible sorting 224 ! i,e. xx,xy,yy,xz,yz,zz 225 !$acc loop vector private(xk1,xl1,xk2,xl2,tii,pqm,pdmk,pdml,kl,kj) 226 do k = 1,3 227 xk1 = ra(k) 228 pdmk = pij*dpint(k,i,i) 229 tii = xk1*ps-pdmk 230 !$acc atomic 231 dipm(k,ii) = dipm(k,ii)+tii 232 ! off-diagonal 233 do l = 1,k-1 234 kl = k*(k-1)/2+l 235 kj = k+l+1 ! the qpint is stored as xx,yy,zz,xy,xz,yz (from integral routine) 236 xl1 = ra(l) 237 pdml = pij*dpint(l,i,i) 238 pqm = pij*qpint(kj,i,i) 239 tii = pdmk*xl1+pdml*xk1-xl1*xk1*ps-pqm 240 !$acc atomic 241 qp(kl,ii) = qp(kl,ii)+tii 242 enddo 243 !diagonal 244 kl = k*(k+1)/2 245 pqm = pij*qpint(k,i,i) 246 tii = 2.0_wp*pdmk*xk1-xk1*xk1*ps-pqm 247 !$acc atomic 248 qp(kl,ii) = qp(kl,ii)+tii 249 enddo 250 enddo 251 !$acc end parallel 252 253 !$acc exit data copyout(dipm(:, :), qp(:, :)) 254 255 ! remove trace 256 do i = 1,nat 257 tii = qp(1,i)+qp(3,i)+qp(6,i) 258 tii = 0.50_wp*tii 259 qp(1:6,i) = 1.50_wp*qp(1:6,i) 260 qp(1,i) = qp(1,i)-tii 261 qp(3,i) = qp(3,i)-tii 262 qp(6,i) = qp(6,i)-tii 263 enddo 264 265 !$acc exit data delete(nao, nat, aoat2(:), s(:, :), p(:, :), dpint(:, :, :), & 266 !$acc& qpint(:, :, :),xyz(:, :)) 267 268end subroutine mmompop 269 270 271! distributed atomic multipole moment interactions: all interactions up to r**-3 272! energy evaluation 273! nat : # of atoms 274! xyz(3,nat) : cartesian coordinates 275! q(nat) : atomic partial charges 276! dipm(3,nat) : cumulative atomic dipole moments (x,y,z) 277! qp(6,nat) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 278! gab3,gab5 : damped R**-3 and R**-5 Coulomb laws, dimension: nat*(nat+1)/2 279! multiplication with numerator then leads to R**-2 and R**-3 decay, respectively 280! e : E_AES 281subroutine aniso_electro(aesData,nat,at,xyz,q,dipm,qp,gab3,gab5,e,epol) 282 use xtb_lin, only : lin 283 implicit none 284 class(TMultipoleData), intent(in) :: aesData 285 integer, intent(in) :: nat,at(:) 286 real(wp), intent(in) :: xyz(:,:),q(:) 287 real(wp), intent(inout) :: e 288 real(wp) qp1(6),rr(3),dp1(3),rij(3) 289 real(wp) edd,e01,e02,e11,r2,tt,tt3,q1,qs2 290 real(wp) ed,eq,epol 291 ! stuff for potential 292 real(wp), intent(in) :: gab3(:,:),gab5(:,:) 293 real(wp), intent(in) :: dipm(:,:),qp(:,:) 294 integer, parameter :: idx(3, 3) = reshape([1, 2, 4, 2, 3, 5, 4, 5, 6], [3, 3]) 295 296 integer i,j,k,l,m,ki,kj,kl 297 298 ! acc enter data copyin(at, xyz, q, dipm, qp, gab3, gab5, & 299 ! acc& aesData, aesData%dipKernel(:), aesData%quadKernel(:)) 300 301 ! acc kernels 302 e = 0.0_wp 303 epol = 0.0_wp 304 e01 = 0.0_wp 305 e02 = 0.0_wp 306 e11 = 0.0_wp 307 ! acc end kernels 308 309 ! acc parallel private(qp1, rr, dp1, rij) 310 311 ! acc loop gang 312 do i = 1, nat 313 q1 = q(i) 314 rr(1:3) = xyz(1:3,i) 315 dp1(1:3) = dipm(1:3,i) 316 qp1(1:6) = qp(1:6,i) 317 ! test: semilocal CT correction 318 ! dipole 319 tt = dp1(1)*dp1(1)+dp1(2)*dp1(2)+dp1(3)*dp1(3) 320 ! qpole 321 tt3 = 0.0_wp 322 ! acc loop seq 323 do k = 1,3 324 ! acc loop seq 325 do l = 1,3 326 kl = idx(l,k) 327 tt3 = tt3+qp1(kl)*qp1(kl) 328 enddo 329 enddo 330 eq = aesData%dipKernel(at(i))*tt+tt3*aesData%quadKernel(at(i)) 331 ! acc atomic 332 epol = epol+eq 333 ! --- 334 enddo 335 336 ! acc loop gang collapse(2) 337 do i = 1, nat 338 do j = 1, nat 339 if (j >= i) cycle 340 q1 = q(i) 341 rr(1:3) = xyz(1:3,i) 342 dp1(1:3) = dipm(1:3,i) 343 qp1(1:6) = qp(1:6,i) 344 kj = i*(i-1)/2 + j 345 rij(1:3) = xyz(1:3,j)-rr(1:3) 346 r2 = sum(rij*rij) 347 ed = 0.0_wp 348 eq = 0.0_wp 349 edd = 0.0_wp 350 ! dipole - charge 351 ! acc loop seq 352 do k = 1,3 353 ed = ed+q(j)*dp1(k)*rij(k) 354 ed = ed-dipm(k,j)*q1*rij(k) 355 ! dip-dip & charge-qpole 356 ! acc loop seq 357 do l = 1,3 358 kl = idx(l,k) 359 tt = rij(l)*rij(k) 360 tt3 = 3.0_wp*tt 361 eq = eq+q(j)*qp1(kl)*tt 362 eq = eq+qp(kl,j)*q1*tt 363 edd = edd-dipm(k,j)*dp1(l)*tt3 364 enddo 365 ! diagonal dip-dip term 366 edd = edd+dipm(k,j)*dp1(k)*r2 367 enddo 368 ! acc atomic 369 e01 = e01+ed*gab3(j,i) 370 ! acc atomic 371 e02 = e02+eq*gab5(j,i) 372 ! acc atomic 373 e11 = e11+edd*gab5(j,i) 374 enddo 375 enddo 376 ! acc end parallel 377 ! acc kernels 378 e = e01 + e02 + e11 379 ! acc end kernels 380 ! write(*,'(''d,q,dd'',3f9.5)') e01,e02,e11 381 ! write(*,*) ' semilocal CT corr.: ',epol 382 ! acc exit data delete(aesData, aesData%dipKernel(:), aesData%quadKernel(:), & 383 ! acc& at, xyz, q, dipm, qp, gab3, gab5) 384 385end subroutine aniso_electro 386 387! aniso-electro from Fock matrix elements 388! nao : # of spherical AOs (SAOs) 389! s(nao,nao) : overlap matrix 390! aoat2(nao) : SAO to atom intex 391! dpint : dipole integral matrix, dimension 3,nao*(nao+1)/2 392! qpint : quadrupole integral matrix, dimension 6,nao*(nao+1)/2 393! p(nao,nao) : density matrix 394! vs(nat) : overlap proportional potential 395! vd(3,nat) : dipint proportional potential 396! vq(6,nat) : quadrupole proportional potential 397subroutine fockelectro(nat,nao,aoat2,p,s,dpint,qpint,vs,vd,vq,e) 398 use xtb_lin, only : lin 399 implicit none 400 integer, intent(in) :: nat,nao,aoat2(nao) 401 real(wp), intent(in) :: dpint(3,nao,nao),s(nao,nao) 402 real(wp), intent(in) :: qpint(6,nao,nao),p(nao,nao) 403 real(wp), intent(in) :: vs(nat),vd(3,nat),vq(6,nat) 404 real(wp), intent(out) :: e 405 real(wp) eaes,pji,fji 406 integer i,j,k,l,ii,jj,ij,kl,kj 407 ! CAMM 408 eaes = 0.0_wp 409 ij = 0 410 do i = 1,nao 411 ii = aoat2(i) 412 do j = 1,nao 413 ij = lin(j,i) 414 jj = aoat2(j) 415 fji = 0.0_wp 416 pji = p(j,i) 417 fji = fji+s(j,i)*(vs(ii)+vs(jj)) 418 do k = 1,3 419 fji = fji+dpint(k,i,j)*(vd(k,ii)+vd(k,jj)) 420 enddo 421 do k = 1,6 422 fji = fji+qpint(k,i,j)*(vq(k,ii)+vq(k,jj)) 423 enddo 424 eaes = eaes+pji*fji 425 enddo 426 enddo 427 eaes = 0.250_wp*eaes 428 ! write(*,*) 'EAES',eaes 429 e = eaes 430end subroutine fockelectro 431 432 433! set-up potential terms v, which are proportional to s, d, or q-integrals 434! it is to be multiplied with Sji when stting up Fji (hence, termed vs) 435! comes essentially at no cost, once cumulative atomic quadrupole moments are available. 436! NEW: the CAMMs were already scaled by scalecamm, but the corresponding potential terms 437! including shift terms need to be scaled 438! nat : # of atoms 439! at(nat) : atom-to-element identifier 440! xyz(3,nat) : cartesian coordinates 441! q(nat) : atomic partial charges 442! dipm(3,nat) : atomic dipole moments 443! qp(6,nat) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 444! gab3,gab5 : damped Coulomb laws, dimension: nat*(nat+1)/2 445! vs(nat) : s-proportional potential from all atoms acting on atom i 446! vd(3,nat) : dipint-proportional potential from all atoms acting on atom i 447! vq(6,nat) : qpole-int proportional potential from all atoms acting on atom i 448subroutine setvsdq(aesData,nat,at,xyz,q,dipm,qp,gab3,gab5,vs,vd,vq) 449 use xtb_lin, only : lin 450 implicit none 451 class(TMultipoleData), intent(in) :: aesData 452 integer, intent(in) :: nat,at(:) 453 real(wp), intent(in) :: q(:),dipm(:,:) 454 real(wp), intent(in) :: xyz(:,:),qp(:,:) 455 real(wp), intent(in) :: gab3(:,:) 456 real(wp), intent(in) :: gab5(:,:) 457 real(wp), intent(out) :: vs(:),vd(:,:),vq(:,:) 458 real(wp) ra(3),dra(3),rb(3),stmp,dum3a,dum5a,t1a,t2a,t3a,t4a,r2a 459 real(wp) r2ab,t1b,t2b,t3b,t4b,dum3b,dum5b,dtmp(3),qtmp(6),g3,g5 460 real(wp) qs1,qs2 461 integer i,j,k,l1,l2,ll,m,mx,ki,kj 462 vs = 0.0_wp 463 vd = 0.0_wp 464 vq = 0.0_wp 465 ! set up overlap proportional potential 466 do i = 1,nat 467 ra(1:3) = xyz(1:3,i) 468 stmp = 0.0_wp 469 dtmp = 0.0_wp 470 qtmp = 0.0_wp 471 do j = 1,nat 472 g3 = gab3(j,i) 473 g5 = gab5(j,i) 474 rb(1:3) = xyz(1:3,j) 475 dra(1:3) = ra(1:3)-rb(1:3) 476 dum3a = 0.0_wp ! collect gab3 dependent terms 477 dum5a = 0.0_wp ! collect gab5 dependent terms 478 r2a = 0.0_wp 479 r2ab = 0.0_wp 480 t1a = 0.0_wp 481 t2a = 0.0_wp 482 t3a = 0.0_wp 483 t4a = 0.0_wp 484 ll = 0 485 do l1 = 1,3 486 ! potential from dipoles 487 r2a = r2a+ra(l1)*ra(l1) ! R_C * R_C 488 r2ab = r2ab+dra(l1)*dra(l1) ! R_AC * R_AC 489 t1a = t1a+ra(l1)*dra(l1) ! R_C * R_AC : for dip-q (q-shift) and dip-dip (q-shift) 490 t2a = t2a+dipm(l1,j)*dra(l1) ! mu_A * R_AC : for q-dip and dip-dip (q-shift) 491 t3a = t3a+ra(l1)*dipm(l1,j) ! R_C * mu_A : for diag. dip-dip (q-shift) 492 t4a = t4a+dra(l1)*dra(l1)*ra(l1)*ra(l1) ! (R_C o R_AC)**"2(square of Hadamard product) : 493 ! results from trace remove from q-pole (q-shift) 494 do l2 = 1,3 495 ll = lin(l1,l2) 496 ! potential from quadrupoles 497 dum5a = dum5a-qp(ll,j)*dra(l1)*dra(l2) & 498 & -1.50_wp*q(j)*dra(l1)*dra(l2)*ra(l1)*ra(l2) 499 if(l2.ge.l1) cycle 500 ki = l1+l2+1 501 qtmp(ki) = qtmp(ki)-3.0_wp*q(j)*g5*dra(l2)*dra(l1) 502 enddo 503 qtmp(l1) = qtmp(l1)-1.50_wp*q(j)*g5*dra(l1)*dra(l1) 504 enddo 505 ! 506 ! set up S-dependent potential 507 dum3a = -t1a*q(j)-t2a ! dip-q (q-shift) and q-dip 508 dum5a = dum5a+t3a*r2ab-3.0_wp*t1a*t2a & !dip-dip (q-shift terms) 509 & +0.50_wp*q(j)*r2a*r2ab !qpole-q (q-shift, trace removal) 510 stmp = stmp+dum5a*g5+dum3a*g3 511 do l1 = 1,3 512 dum3a = dra(l1)*q(j) 513 dum5a = 3.0_wp*dra(l1)*t2a & ! dipint-dip 514 & -r2ab*dipm(l1,j) & ! dipint-dip (diagonal) 515 & -q(j)*r2ab*ra(l1) & ! qpole-q (dipint-shift, trace removal) 516 & +3.0_wp*q(j)*dra(l1)*t1a ! qpole-q (dipint-shift) 517 dtmp(l1) = dtmp(l1)+dum3a*g3+dum5a*g5 518 qtmp(l1) = qtmp(l1)+0.50_wp*r2ab*q(j)*g5 !remove trace term 519 enddo 520 enddo 521 vs(i) = stmp ! q terms 522 vd(1:3,i) = dtmp(1:3) ! dipints from atom i 523 vq(1:6,i) = qtmp(1:6) ! qpints from atom i 524 ! --- CT correction terms 525 qs1 = aesData%dipKernel(at(i))*2.0_wp 526 qs2 = aesData%quadKernel(at(i))*6.0_wp ! qpole pot prefactors 527 t3a = 0.0_wp 528 t2a = 0.0_wp 529 do l1 = 1,3 530 ! potential from dipoles 531 t3a = t3a+ra(l1)*dipm(l1,i)*qs1 ! R_C * mu_C : for diag. dip-dip 532 vd(l1,i) = vd(l1,i)-qs1*dipm(l1,i) 533 do l2 = 1,l1-1 534 ! potential from quadrupoles 535 ll = lin(l1,l2) 536 ki = l1+l2+1 537 vq(ki,i) = vq(ki,i)-qp(ll,i)*qs2 538 t3a = t3a-ra(l1)*ra(l2)*qp(ll,i)*qs2 539 vd(l1,i) = vd(l1,i)+ra(l2)*qp(ll,i)*qs2 540 vd(l2,i) = vd(l2,i)+ra(l1)*qp(ll,i)*qs2 541 enddo 542 ! diagonal 543 ll = lin(l1,l1) 544 vq(l1,i) = vq(l1,i)-qp(ll,i)*qs2*0.50_wp 545 t3a = t3a-ra(l1)*ra(l1)*qp(ll,i)*qs2*0.50_wp 546 vd(l1,i) = vd(l1,i)+ra(l1)*qp(ll,i)*qs2 547 ! collect trace removal terms 548 t2a = t2a+qp(ll,i) 549 enddo 550 vs(i) = vs(i)+t3a 551 ! trace removal 552 t2a = t2a*aesData%quadKernel(at(i)) 553 do l1 = 1,3 554 vq(l1,i) = vq(l1,i)+t2a 555 vd(l1,i) = vd(l1,i)-2.0_wp*ra(l1)*t2a 556 vs(i) = vs(i)+t2a*ra(l1)*ra(l1) 557 enddo 558 ! --- 559 enddo 560 561 ! call prmat(6,vs,nat,0,'vs') 562 563end subroutine setvsdq 564 565! set-up potential terms v used for nuclear gradients 566! here, the shift terms are removed, since we use multipole derivatives w/ origin at the atoms 567! nat : # of atoms 568! xyz(3,nat) : cartesian coordinates 569! q(nat) : atomic partial charges 570! dipm(3,nat) : atomic dipole moments 571! qp(6,nat) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 572! gab3,gab5 : damped Coulomb laws, dimension: nat*(nat+1)/2 573! vs(nat) : s-proportional potential from all atoms acting on atom i 574! vd(3,nat) : dipint-proportional potential from all atoms acting on atom i 575! vq(6,nat) : qpole-int proportional potential from all atoms acting on atom i 576subroutine setdvsdq(aesData,nat,at,xyz,q,dipm,qp,gab3,gab5,vs,vd,vq) 577 use xtb_lin, only : lin 578 implicit none 579 class(TMultipoleData), intent(in) :: aesData 580 integer, intent(in) :: nat,at(:) 581 real(wp), intent(in) :: q(:),dipm(:,:) 582 real(wp), intent(in) :: xyz(:,:),qp(:,:) 583 real(wp), intent(in) :: gab3(:,:) 584 real(wp), intent(in) :: gab5(:,:) 585 real(wp), intent(out) :: vs(:),vd(:,:),vq(:,:) 586 real(wp) ra(3),dra(3),rb(3),stmp,dum3a,dum5a,t1a,t2a,t3a,t4a,r2a 587 real(wp) r2ab,t1b,t2b,t3b,t4b,dum3b,dum5b,dtmp(3),qtmp(6),g3,g5 588 real(wp) qs1,qs2 589 integer i,j,k,l1,l2,ll,m,mx,ki,kj 590 vs = 0.0_wp 591 vd = 0.0_wp 592 vq = 0.0_wp 593 ! set up overlap proportional potential 594 do i = 1,nat 595 ra(1:3) = xyz(1:3,i) 596 stmp = 0.0_wp 597 dtmp = 0.0_wp 598 qtmp = 0.0_wp 599 do j = 1,nat 600 g3 = gab3(j,i) 601 g5 = gab5(j,i) 602 rb(1:3) = xyz(1:3,j) 603 dra(1:3) = ra(1:3)-rb(1:3) 604 dum3a = 0.0_wp ! collect gab3 dependent terms 605 dum5a = 0.0_wp ! collect gab5 dependent terms 606 r2a = 0.0_wp 607 r2ab = 0.0_wp 608 t2a = 0.0_wp 609 ll = 0 610 do l1 = 1,3 611 ! potential from dipoles 612 r2ab = r2ab+dra(l1)*dra(l1) ! R_AC * R_AC 613 t2a = t2a+dipm(l1,j)*dra(l1) ! mu_A * R_AC : for q-dip and dip-dip (q-shift) 614 do l2 = 1,3 615 ll = lin(l1,l2) 616 ! potential from quadrupoles 617 dum5a = dum5a-qp(ll,j)*dra(l1)*dra(l2) 618 if(l2.ge.l1) cycle 619 ki = l1+l2+1 620 qtmp(ki) = qtmp(ki)-3.0_wp*q(j)*g5*dra(l2)*dra(l1) 621 enddo 622 qtmp(l1) = qtmp(l1)-1.50_wp*q(j)*g5*dra(l1)*dra(l1) 623 enddo 624 dum3a = -t2a ! q-dip ! w/o shift terms 625 stmp = stmp+dum3a*g3+dum5a*g5 626 do l1 = 1,3 627 dum3a = dra(l1)*q(j) ! w/o shift terms 628 dum5a = 3.0_wp*dra(l1)*t2a & ! dipint-dip 629 & -r2ab*dipm(l1,j) ! dipint-dip (diagonal) 630 dtmp(l1) = dtmp(l1)+dum3a*g3+dum5a*g5 631 qtmp(l1) = qtmp(l1)+0.50_wp*r2ab*q(j)*g5 !remove trace term 632 enddo 633 enddo 634 vs(i) = stmp 635 vd(1:3,i) = dtmp(1:3) 636 vq(1:6,i) = qtmp(1:6) 637 ! --- CT correction terms 638 qs1 = aesData%dipKernel(at(i))*2.0_wp 639 qs2 = aesData%quadKernel(at(i))*6.0_wp ! qpole pot prefactors 640 t2a = 0.0_wp 641 do l1 = 1,3 642 ! potential from dipoles 643 vd(l1,i) = vd(l1,i)-qs1*dipm(l1,i) 644 do l2 = 1,l1-1 645 ! potential from quadrupoles 646 ll = lin(l1,l2) 647 ki = l1+l2+1 648 vq(ki,i) = vq(ki,i)-qp(ll,i)*qs2 649 enddo 650 ll = lin(l1,l1) 651 vq(l1,i) = vq(l1,i)-qp(ll,i)*qs2*0.50_wp 652 ! collect trace removal terms 653 t2a = t2a+qp(ll,i) 654 enddo 655 ! trace removal 656 t2a = t2a*aesData%quadKernel(at(i)) 657 do l1 = 1,3 658 vq(l1,i) = vq(l1,i)+t2a 659 enddo 660 ! --- 661 662 663 enddo 664 665 ! call prmat(6,vs,nat,0,'vs') 666 667end subroutine setdvsdq 668 669! molmom: computes molecular multipole moments from CAMM 670! n : # of atoms 671! xyz(3,n) : cartesian coordinates 672! q(n) : atomic partial charges 673! dipm(3,n) : cumulative atomic dipole moments (x,y,z) 674! qp(6,n) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 675subroutine molmom(iunit,n,xyz,q,dipm,qp,dip,d3) 676 use xtb_mctc_convert 677 use xtb_lin, only : lin 678 implicit none 679 integer, intent(in) :: iunit 680 integer, intent(in) :: n 681 real(wp), intent(in) :: xyz(:,:),q(:),dipm(:,:),qp(:,:) 682 real(wp), intent(out) :: dip,d3(:) 683 real(wp) rr1(3),rr2(3),tma(6),tmb(6),tmc(6),dum 684 integer i,j,k,l 685 rr1 = 0.0_wp 686 rr2 = 0.0_wp 687 write(iunit,'(a)') 688 do i = 1,n 689 do j = 1,3 690 rr1(j) = rr1(j)+q(i)*xyz(j,i) 691 rr2(j) = rr2(j)+dipm(j,i) 692 enddo 693 enddo 694 d3(1:3) = rr1(1:3)+rr2(1:3) 695 dip = sqrt(d3(1)**2+d3(2)**2+d3(3)**2) 696 write(iunit,'(a)',advance='yes')'molecular dipole:' 697 write(iunit,'(a)',advance='no')' ' 698 write(iunit,'(a)',advance='yes') & 699 & 'x y z tot (Debye)' 700 write(iunit,'(a,3f12.3)') ' q only: ',rr1(1:3) 701 write(iunit,'(a,4f12.3)') ' full: ',d3(1:3),dip*autod 702 703 tma = 0.0_wp 704 tmb = 0.0_wp 705 tmc = 0.0_wp 706 do i = 1,n 707 l = 0 708 do j = 1,3 709 do k = 1,j 710 l = lin(k,j) 711 tma(l) = tma(l)+xyz(j,i)*xyz(k,i)*q(i) 712 tmb(l) = tmb(l)+dipm(k,i)*xyz(j,i)+dipm(j,i)*xyz(k,i) 713 tmc(l) = tmc(l)+qp(l,i) 714 enddo 715 enddo 716 enddo 717 ! remove traces and multiply with 3/2 in q and dip parts 718 dum = tma(1)+tma(3)+tma(6) 719 dum = 0.50_wp*dum 720 tma = 1.50_wp*tma 721 l = 0 722 do j = 1,3 723 l = l+j 724 tma(l) = tma(l)-dum 725 enddo 726 dum = tmb(1)+tmb(3)+tmb(6) 727 dum = 0.50_wp*dum 728 tmb = 1.50_wp*tmb 729 l = 0 730 do j = 1,3 731 l = l+j 732 tmb(l) = tmb(l)-dum 733 enddo 734 write(iunit,'(a)',advance='yes')'molecular quadrupole (traceless):' 735 write(iunit,'(a)',advance='no')' ' 736 write(iunit,'(a)',advance='no')'xx xy yy ' 737 write(iunit,'(a)',advance='yes')'xz yz zz' 738 write(iunit,'(a,6f12.3)') ' q only: ',tma(1:6) 739 write(iunit,'(a,6f12.3)') ' q+dip: ',tma(1:6)+tmb(1:6) 740 write(iunit,'(a,6f12.3)') ' full: ',tma(1:6)+tmb(1:6)+tmc(1:6) 741 742end subroutine molmom 743 744! gradient evaluation from 745! cumulative atomic multipole moment interactions: all interactions up to r**-3 746! nat : # of atoms 747! xyz(3,nat) : cartesian coordinates 748! q(nat) : atomic partial charges 749! dipm(3,nat) : cumulative atomic dipole moments (x,y,z) 750! qp(6,nat) : traceless(!) cumulative atomic quadrupole moments (xx,xy,yy,xz,yz,zz) 751! gab3,gab5 : damped R**-3 and R**-5 Coulomb laws, dimension: nat*(nat+1)/2 752! multiplication with numerator then leads to R**-2 and R**-3 decay, respectively 753! radcn(nat) : CN-depentent atomic radii 754! dcn(3,i,j) : derivative of radcn(j) w.r.t. cartesian directions of i 755! exj : exponent in gab, gab3, and gab5 - determines interpolation 756! g : nuclear gradient (3) 757subroutine aniso_grad(nat,at,xyz,q,dipm,qp,kdmp3,kdmp5, & 758 & radcn,dcn,gab3,gab5,g) 759 use xtb_lin, only : lin 760 !gab3 Hellmann-Feynman terms correct, shift terms to be tested yet 761 implicit none 762 integer, intent(in) :: nat,at(:) 763 real(wp), intent(in) :: xyz(:,:),q(:),dipm(:,:),qp(:,:) 764 real(wp), intent(in) :: gab3(:,:),gab5(:,:) 765 real(wp), intent(in) :: kdmp3,kdmp5,radcn(:),dcn(:,:,:) 766 real(wp), intent(inout) :: g(:,:) 767 real(wp) qp1(6),rr(3),dip(3),rij(3) 768 real(wp) ed,eq,edd,e01,e02,e11,r2,tt,tt3,q1,dxi 769 real(wp) tmp2,tmp3,rab,rabi,ddm2,ddm3a,ddm3b,qqa,qqb 770 real(wp) dgab3,dgab5,damp1,damp2,ddamp,qs2 771 772 integer i,j,k,l,m,ki,kj,kl 773 do i = 1,nat 774 q1 = q(i) 775 rr(1:3) = xyz(1:3,i) 776 dip(1:3) = dipm(1:3,i) 777 qp1(1:6) = qp(1:6,i) 778 tmp2 = 0.0_wp ! cumulate terms propto CN gradient - to scale only quadratically 779 do j = 1,nat ! loop over other atoms 780 if(i.eq.j) cycle 781 kj = lin(j,i) 782 rij(1:3) = xyz(1:3,j)-rr(1:3) 783 r2 = sum(rij*rij) 784 rabi = 1.0_wp/sqrt(r2) 785 ! call dzero(2.0_wp,rabi,at(i),at(j),damp,ddamp) 786 call dzero(kdmp3,rabi,radcn(i),radcn(j),damp1,ddamp) 787 dgab3 = dgab(3.0_wp,rabi,damp1,ddamp) 788 ! call dzero(3.0_wp,rabi,at(i),at(j),damp,ddamp) 789 call dzero(kdmp5,rabi,radcn(i),radcn(j),damp2,ddamp) 790 dgab5 = dgab(5.0_wp,rabi,damp2,ddamp) 791 !!! DEBUG 792 ! dgab3 = 0.0_wp 793 ! dgab5 = 0.0_wp 794 ! dgab3 = dgab3*100.0_wp 795 ! dgab5 = dgab5*100.0_wp 796 !!! 797 ed = 0.0_wp 798 edd = 0.0_wp 799 eq = 0.0_wp 800 ! dipole - charge 801 do k = 1,3 802 ed = ed+q(j)*dip(k)*rij(k) 803 ed = ed-dipm(k,j)*q1*rij(k) 804 tt = q1*dipm(k,j)-q(j)*dip(k) 805 ! part of dip-q derivative 806 g(k,i) = g(k,i)+gab3(j,i)*tt 807 ! dip-dip & charge-qpole 808 ddm2 = 0.0_wp 809 ddm3a = 0.0_wp 810 ddm3b = 0.0_wp 811 qqa = 0.0_wp 812 qqb = 0.0_wp 813 do l = 1,3 814 kl = lin(l,k) 815 tt = rij(l)*rij(k) 816 tt3 = 3.0_wp*tt 817 eq = eq+q(j)*qp1(kl)*tt 818 eq = eq+qp(kl,j)*q1*tt 819 edd = edd-dipm(k,j)*dip(l)*tt3 820 ! extra d-d terms 821 ddm2 = ddm2+dipm(l,j)*dip(l) 822 ddm3a = ddm3a+dip(l)*rij(l) 823 ddm3b = ddm3b+dipm(l,j)*rij(l) 824 ! extra q-qpole terms 825 qqa = qqa+rij(l)*qp(kl,j) 826 qqb = qqb+rij(l)*qp1(kl) 827 enddo 828 edd = edd+dipm(k,j)*dip(k)*r2 829 g(k,i) = g(k,i)-2.0_wp*gab5(j,i)*ddm2*rij(k) 830 g(k,i) = g(k,i)+3.0_wp*gab5(j,i)*ddm3a*dipm(k,j) 831 g(k,i) = g(k,i)+3.0_wp*gab5(j,i)*ddm3b*dip(k) 832 g(k,i) = g(k,i)-2.0_wp*gab5(j,i)*qqa*q1 833 g(k,i) = g(k,i)-2.0_wp*gab5(j,i)*qqb*q(j) 834 enddo 835 do k = 1,3 836 dxi = rij(k)*rabi 837 g(k,i) = g(k,i)-ed*dgab3*dxi 838 g(k,i) = g(k,i)-(eq+edd)*dgab5*dxi 839 enddo 840 ! collect terms for CN-dependent part 841 rab = 0.50_wp*(radcn(i)+radcn(j)) 842 tmp3 = ed*kdmp3*gab3(j,i)*(damp1/rab)*(rab*rabi)**kdmp3 843 tmp2 = tmp2+tmp3 844 tmp3 = (eq+edd)*kdmp5*gab5(j,i)*(damp2/rab)*(rab*rabi)**kdmp5 845 tmp2 = tmp2+tmp3 846 enddo 847 ! CN-dependent part - O(N^2) 848 tmp2 = 3.0_wp*tmp2 849 g(:,:) = g-tmp2*dcn(:,:,i) 850 851 enddo 852end subroutine aniso_grad 853 854 855! check and print sparsity w.r.t. individual contribution 856! to get an idea 857subroutine checkspars(nao,ndp,nqp,nmat,matlist,mqlst,mdlst) 858 use xtb_lin, only : lin 859 implicit none 860 integer,intent(in) :: ndp,nqp,nmat,nao 861 integer,intent(in) :: matlist(:,:) 862 integer,intent(in) :: mqlst(:,:),mdlst(:,:) 863 integer :: i,j,m,k,ntot,mi,mj,ki,kj,mm,kk 864 logical,allocatable :: nzero(:) 865 ! check overall sparsity 866 allocate(nzero(nao*(nao+1)/2)) 867 nzero = .false. 868 do k = 1,ndp 869 ki = mdlst(1,k) 870 kj = mdlst(2,k) 871 kk = lin(ki,kj) 872 nzero(kk) = .true. 873 enddo 874 do k = 1,nqp 875 ki = mqlst(1,k) 876 kj = mqlst(2,k) 877 kk = lin(ki,kj) 878 nzero(kk) = .true. 879 enddo 880 do k = 1,nmat 881 ki = matlist(1,k) 882 kj = matlist(2,k) 883 kk = lin(ki,kj) 884 nzero(kk) = .true. 885 enddo 886 mm = nao*(nao+1)/2 887 ntot = 0 888 do i = 1,mm 889 if(nzero(i)) ntot = ntot+1 890 enddo 891 write(*,'(a)',advance='yes') ' ' 892 write(*,'(a)')'% of non-zero elements in H:' 893 write(*,'('' by overlap:'',f6.2)') & 894 & 100.*float(nmat)/float(mm) 895 write(*,'('' by dipole ints.:'',f6.2)') & 896 & 100.*float(ndp)/float(mm) 897 write(*,'('' by quadrupole ints.:'',f6.2)') & 898 & 100.*float(nqp)/float(mm) 899 write(*,'('' total:'',f6.2)') & 900 & 100.*float(ntot)/float(mm) 901 write(*,'(a)',advance='yes') ' ' 902 deallocate(nzero) 903end subroutine checkspars 904 905! zero-damped gab 906subroutine mmomgabzero(nat,at,xyz,kdmp3,kdmp5,radcn,gab3,gab5) 907 implicit none 908 integer, intent(in) :: nat,at(:) 909 real(wp), intent(in) :: xyz(:,:),radcn(:) 910 real(wp), intent(in) :: kdmp3,kdmp5 911 real(wp), intent(out) :: gab3(:,:),gab5(:,:) 912 real(wp) damp,ddamp 913 914 real(wp) tmp1,tmp2,rr(3) 915 integer i,j,k,l,lin 916 917 !!!!!!! set up damped Coulomb operators for multipole interactions 918 gab3 = 0.0_wp ! for r**-2 decaying q-dip term 919 gab5 = 0.0_wp ! for r**-3 decaying terms (q-qpol,dip-dip) 920 do i = 1,nat 921 do j = 1,nat 922 if (j.ge.i) cycle 923 rr(1:3) = xyz(1:3,i) 924 l = i*(i-1)/2+j 925 tmp2 = 0.0_wp 926 do k = 1,3 927 tmp1 = xyz(k,j)-rr(k) 928 tmp2 = tmp2+tmp1**2 929 enddo 930 tmp1 = 1.0_wp/sqrt(tmp2) 931 ! call dzero(2.0_wp,tmp1,at(i),at(j),damp,ddamp) 932 call dzero(kdmp3,tmp1,radcn(i),radcn(j),damp,ddamp) 933 gab3(j,i) = gab(3.0_wp,tmp1,damp) 934 gab3(i,j) = gab3(j,i) 935 ! call dzero(3.0_wp,tmp1,at(i),at(j),damp,ddamp) 936 call dzero(kdmp5,tmp1,radcn(i),radcn(j),damp,ddamp) 937 gab5(j,i) = gab(5.0_wp,tmp1,damp) 938 gab5(i,j) = gab5(j,i) 939 enddo 940 enddo 941 942end subroutine mmomgabzero 943 944 945!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 946! setup CN-dependet atomic radii 947! n : # of atoms 948! at(n) : atomic number array 949! cn(n) : coordination number of atoms 950! shift : global offset from cnval (parameter) 951! expo : exponent scaling/ steepness (parameter) 952! rmax : maximum radius (parameter) 953! radcn : CN-dependent radius 954subroutine get_radcn(aesData,n,at,cn,shift,expo,rmax,radcn) 955 implicit none 956 class(TMultipoleData), intent(in) :: aesData 957 integer, intent (in) :: n,at(:) 958 real(wp), intent (in) :: cn(:),shift,expo,rmax 959 real(wp), intent (out) :: radcn(:) 960 real(wp) rco,t1,t2 961 integer i,j 962 do i = 1,n 963 rco = aesData%multiRad(at(i)) ! base radius of element 964 t1 =cn(i)-aesData%valenceCN(at(i))-shift ! CN - VALCN - SHIFT 965 t2 =rco +(rmax-rco)/(1.0_wp+exp(-expo*t1)) 966 radcn(i) = t2 967 enddo 968end subroutine get_radcn 969!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 970! derivative of CN-dependet atomic radii w.r.t. other atoms 971! n : # of atoms 972! at(n) : atomic number array 973! cn(n) : coordination number of atoms 974! shift : global offset from cnval (parameter) 975! expo : exponent scaling/ steepness (parameter) 976! rmax : maximum radius (parameter) 977! dcn : on input : derivatives of CN(j) w.r.t. Cart. directions of i 978! : on output : derivatives of RADCN(j) w.r.t. Cart. directions of i 979subroutine dradcn(aesData,n,at,cn,shift,expo,rmax,dcn) 980 implicit none 981 class(TMultipoleData), intent(in) :: aesData 982 integer, intent (in) :: n,at(:) 983 real(wp), intent (in) :: cn(:),shift,expo,rmax 984 real(wp), intent (inout) :: dcn(:,:,:) 985 real(wp) rco,t1,t2,t3,t4,tmp1,tmp2 986 integer i,j,k 987 do i = 1,n 988 rco = aesData%multiRad(at(i)) ! base radius of element 989 t1 =exp(-expo*(cn(i)-aesData%valenceCN(at(i))-shift)) ! CN - VALCN - SHIFT 990 t2 =(rmax-rco)/(1.0_wp+2.0_wp*t1+t1*t1) 991 t2 = t2*expo*t1 992 dcn(:,:,i) = dcn(:,:,i)*t2 993 end do 994end subroutine dradcn 995 996!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 997! zero-damping function and derivative 998! rscal : scaling of radii 999! dex : exponent in zero-damping function 1000! rabinv : inverse distance, i.e., 1/Rab 1001! aradi,radj : CN-dependent atomic number radii 1002! damp : zero damping function 1003! ddamp : derivative w.r.t. Rab 1004subroutine dzero(dex,rabinv,radi,radj,damp,ddamp) 1005 implicit none 1006 real(wp), intent (in) :: radi,radj 1007 real(wp), intent (in) :: rabinv,dex 1008 real(wp), intent (out) :: damp,ddamp 1009 real(wp) rco,f1,f2 1010 ! rco = 2.0/(1./aesData%multiRad(ati)+1./aesData%multiRad(atj)) ! unstable 1011 ! rco = sqrt(aesData%multiRad(ati)*aesData%multiRad(atj)) ! unstable 1012 rco = 0.5*(radi+radj) 1013 ! zero-damping function and gradient w.r.t. Rab 1014 damp = 1.0_wp/(1.0_wp+6.0_wp*(rco*rabinv)**dex) 1015 ddamp = -dex*rabinv*(damp*damp-damp) 1016end subroutine dzero 1017 1018!!! gab - computes the damped Coulomb type interaction with dex decay: 1019! gab = damp * Rab**(-dex) 1020! dex : exponent defining the decay: Rab**(-dex) 1021! rabinv : inverse distance, i.e., 1/Rab 1022! damp : zero damping function 1023real(wp) function gab(dex,rabinv,damp) 1024 implicit none 1025 real(wp) dex,rabinv,damp 1026 ! compute r**dex decaying intermediate 1027 gab = damp*(rabinv**dex) ! LR-decay * damping 1028end function gab 1029 1030 1031!!! dgab - computes the derivative w.r.t. Rab of a damped Coulomb type interaction with dex decay: 1032! gab = damp * Rab**(-dex) 1033! dex : exponent defining the decay: Rab**(-dex) 1034! rabinv : inverse distance, i.e., 1/Rab 1035! damp : zero damping function 1036! ddamp : derivative of damping function 1037real(wp) function dgab(dex,rabinv,damp,ddamp) 1038 implicit none 1039 real(wp) dex,rabinv,damp,ddamp,tmp1 1040 ! compute r**dex decaying intermediate 1041 tmp1 = -dex*rabinv*(rabinv**dex) ! LR-decay derivative 1042 dgab = tmp1*damp+ddamp*rabinv**dex 1043end function dgab 1044 1045!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1046 1047subroutine gfn2broyden_diff(n,istart,nbr,dipm,qp,q_in,dq) 1048 implicit none 1049 integer, intent (in) :: n,nbr 1050 integer, intent (inout) :: istart 1051 real(wp), intent (in) :: dipm(:,:),qp(:,:),q_in(:) 1052 real(wp), intent (inout) :: dq(:) 1053 integer i,j,k 1054 k = istart 1055 do i = 1,n 1056 do j = 1,3 1057 k = k+1 1058 dq(k) = dipm(j,i)-q_in(k) 1059 enddo 1060 do j = 1,6 1061 k = k+1 1062 dq(k) = qp(j,i)-q_in(k) 1063 enddo 1064 enddo 1065 istart = k 1066 1067end subroutine gfn2broyden_diff 1068 1069subroutine gfn2broyden_save(n,istart,nbr,dipm,qp,q_in) 1070 implicit none 1071 integer, intent (in) :: n,nbr 1072 integer, intent (inout) :: istart 1073 real(wp), intent (in) :: dipm(:,:),qp(:,:) 1074 real(wp), intent (inout) :: q_in(:) 1075 integer i,j,k 1076 k = istart 1077 do i = 1,n 1078 do j = 1,3 1079 k = k+1 1080 q_in(k) = dipm(j,i) 1081 enddo 1082 do j = 1,6 1083 k = k+1 1084 q_in(k) = qp(j,i) 1085 enddo 1086 enddo 1087 istart = k 1088 1089end subroutine gfn2broyden_save 1090 1091 1092 1093subroutine gfn2broyden_out(n,istart,nbr,q_in,dipm,qp) 1094 implicit none 1095 integer, intent (in) :: n,nbr 1096 integer, intent (inout) :: istart 1097 real(wp), intent (in) :: q_in(:) 1098 real(wp), intent (out) :: dipm(:,:),qp(:,:) 1099 integer i,j,k 1100 k = istart 1101 do i = 1,n 1102 do j = 1,3 1103 k = k+1 1104 dipm(j,i) = q_in(k) 1105 enddo 1106 do j = 1,6 1107 k = k+1 1108 qp(j,i) = q_in(k) 1109 enddo 1110 enddo 1111 istart = k 1112end subroutine gfn2broyden_out 1113 1114 1115end module xtb_aespot 1116