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_qmdff 19 use, intrinsic :: iso_fortran_env, only : output_unit 20 use xtb_mctc_accuracy, only : wp 21 implicit none 22! integer, private :: ndim 23! parameter(ndim=50000) 24 integer, private,parameter :: max_elem = 94 25 integer, private,parameter :: ntterm = 4 26 integer, private :: nbond = 0 27 integer, private :: nbond12 = 0 28 integer, private :: nangl = 0 29 integer, private :: ntors = 0 30 integer, private :: nhb = 0 31 integer, private :: nnci = 0 32 integer, private,allocatable :: bond (:,:) 33 integer, private,allocatable :: angl (:,:) 34 integer, private,allocatable :: tors (:,:) 35 integer, private,allocatable :: hb (:,:) 36 integer, private,allocatable :: nci (:,:) 37 real(wp),private,allocatable :: vbond(:,:) 38 real(wp),private,allocatable :: vangl(:,:) 39 real(wp),private,allocatable :: vtors(:,:) 40 real(wp),private,allocatable :: vhb (:,:) 41 real(wp),private :: scalehb(max_elem) = 0.0_wp 42 real(wp),private :: scalexb(max_elem) = 0.0_wp 43 real(wp),private :: morsethr = 99.0_wp 44 real(wp),private,allocatable :: qff(:) 45 real(wp),private,allocatable :: c6ff(:,:) 46 real(wp),private :: eps1(6)=(/0.00_wp,0.00_wp,0.85_wp,1.00_wp,1.00_wp,0.00_wp/) 47 real(wp),private :: eps2(6)=(/0.00_wp,0.00_wp,0.50_wp,0.50_wp,1.00_wp,1.00_wp/) 48 real(wp),private :: zabff (max_elem,max_elem) = 0.0_wp 49 real(wp),private :: r094ff(max_elem,max_elem) = 0.0_wp 50 real(wp),private :: sr42ff(max_elem,max_elem) = 0.0_wp 51 real(wp),private :: r0abff(max_elem,max_elem) = 0.0_wp 52 53 real(wp),private,parameter :: pi = 3.14159265358979323846264338327950_wp 54 real(wp),private,parameter :: pi2 = 6.28318530717958623199592693708837_wp 55 real(wp),private,parameter :: pi6 = 961.389193575304212170602974626298_wp 56 real(wp),private,parameter :: spi = 1.77245385090551599275151910313925_wp 57contains 58 59subroutine ff_ini(n,at,xyz,cn,s6) 60 use xtb_mctc_param, only : r2r4 => sqrt_z_r4_over_r2, & 61 & rcov => covalent_radius_d3 62 implicit none 63 integer, intent(in) :: n 64 integer, intent(in) :: at(n) 65 real(wp),intent(in) :: xyz(3,n) 66 real(wp),intent(in) :: cn(n) 67 real(wp),intent(in) :: s6 68 69 real(wp) :: e 70 real(wp) :: g(3,n) 71 real(wp) :: vz(max_elem) 72 real(wp) :: s8 = 2.70_wp 73 real(wp) :: a1 = 0.45_wp 74 real(wp) :: a2 = 4.00_wp 75 real(wp) :: c6 76 integer :: i,j,i1,i2,iz1,iz2 77 logical :: exist 78 79! include 'ffparam.f90' 80 81 ! clear heap from previous runs 82 if (allocated(bond)) deallocate(bond) 83 if (allocated(angl)) deallocate(angl) 84 if (allocated(tors)) deallocate(tors) 85 if (allocated(hb)) deallocate(hb) 86 if (allocated(nci)) deallocate(nci) 87 if (allocated(vbond)) deallocate(vbond) 88 if (allocated(vangl)) deallocate(vangl) 89 if (allocated(vtors)) deallocate(vtors) 90 if (allocated(vhb)) deallocate(vhb) 91 if (allocated(qff)) deallocate(qff) 92 if (allocated(c6ff)) deallocate(c6ff) 93 94 inquire(file='solvent',exist=exist) 95 if(.not.exist) then 96 call raise('E','FF run requested but solvent file does not exist!') 97 endif 98 99 write(output_unit,'("initializing FF ..")') 100 101! if(n.gt.5000) call raise('E','too many atoms for FF',1) ! fixed this for you 102 103 scalehb=0 104 scalehb(7 )=0.8 105 scalehb(8 )=0.3 106 scalehb(9 )=0.1 107 scalehb(15 )=2.0 108 scalehb(16 )=2.0 109 scalehb(17 )=2.0 110 scalehb(34 )=2.0 111 scalehb(35 )=2.0 112 scalexb=0 113 scalexb(17 )=0.30 114 scalexb(35 )=0.60 115 scalexb(53 )=0.80 116 scalexb(85 )=1.00 117 eps1(1)=0 118 eps1(2)=0 119 eps1(3)=0.85 120 eps1(4)=1 121 eps1(5)=1 122 eps2(1)=0 123 eps2(2)=0 124 eps2(3)=0.5 125 eps2(4)=0.5 126 eps2(5)=1 127 eps2(6)=1.00 128 eps1(6)=0 129 130 call rdsolvff(n,'solvent') 131 132 call valelff(vz) 133 ! D3 134 s8 = s8 / s6 135 136 do i=1,max_elem 137 do j=1,max_elem 138 sr42ff(j,i)=3.0d0*s8*r2r4(i)*r2r4(j) 139 r094ff(j,i)=a1*sqrt(3.0d0*r2r4(i)*r2r4(j))+a2 140 zabff (j,i)=vz(i)*vz(j) 141 enddo 142 enddo 143 call setr0ab(max_elem,0.52917726d0,r0abff) 144 r0abff = 16.5 / r0abff**1.5 145 146 allocate( c6ff(n,n), source = 0.0_wp ) 147 do i1=1,n 148 iz1=at(i1) 149 do i2=1,i1 150 iz2=at(i2) 151 call getc6(iz1,iz2,cn(i1),cn(i2),c6) 152 c6ff(i2,i1)=c6 * s6 ! simulates solvent for s6 < 1 153 c6ff(i1,i2)=c6 * s6 154 enddo 155 enddo 156 157 158! call ff_eg (n,at,xyz,e,g) 159! call ff_nonb(n,at,xyz,e,g) 160! call ff_hb (n,at,xyz,e,g) 161! write(*,*) 'QMDFF energy on input coordinates: ',e 162! write(*,*) 'QMDFF grad on input coordinates: ',sqrt(sum(g**2)) 163! stop 164 165end subroutine ff_ini 166 167!cccccccccccccccccccccccccccccccccccccccccccccccccccc 168! read solvent coord and FF 169!cccccccccccccccccccccccccccccccccccccccccccccccccccc 170 171subroutine rdsolvff(nat,fname) 172 implicit none 173 character(len=*),intent(in) :: fname 174 175 integer, intent(in) :: nat 176 integer :: at(nat) 177 real(wp) :: xyz(3,nat) 178 179 integer :: i,j,n,nn,nterm,idum,imass,ich 180 real(wp) :: qdum,r,xx(10),dum1,dum2 181 character(len=80) atmp 182 183! include 'ffparam.f90' 184 185 open(newunit=ich,file=fname) 186 187 read (ich,'(a)') atmp 188 call readl(atmp,xx,nn) 189 n=idint(xx(1)) 190 if(n.ne.nat) then 191 write(*,*) n,nat 192 stop 'ff read error' 193 endif 194 qdum=xx(2) 195 read (ich,'(a)') atmp 196 write(output_unit,'(a)') 197 write(output_unit,'(a)')'reading QMDFF intramolecular FF:' 198 write(output_unit,'(a)')'method used for charges '//trim(atmp) 199 if(index(atmp,'morsethr=').ne.0) then 200 call readl(atmp,xx,nn) 201 morsethr=xx(nn) 202 write(output_unit,'(a,f12.8)') 'taking morsethr from file ',morsethr 203 else 204 morsethr=99. 205 endif 206 207 allocate(qff(n), source = 0.0_wp) 208 do i=1,n 209 read (ich,*)at(i),xyz(1:3,i),qff(i),imass 210 enddo 211 read (ich,*) nbond,nangl,ntors,nhb,nnci,nbond12 212! if(nbond.gt.ndim.or.nangl.gt.ndim.or.ntors.gt.ndim) & 213! & call raise('E','too many FF terms',1) 214 allocate ( bond(2,nbond), source = 0 ) 215 allocate ( vbond(3,nbond), source = 0.0_wp ) 216 do i=1,nbond 217 read (ich,*)bond(1:2,i),vbond(1:3,i) 218 enddo 219 220 allocate ( angl(3,nangl), source = 0 ) 221 allocate ( vangl(2,nangl), source = 0.0_wp ) 222 do i=1,nangl 223 read (ich,*) angl(1,i),angl(2,i),angl(3,i),vangl(1:2,i) 224 enddo 225 226 allocate ( tors(6,ntors), source = 0 ) 227 allocate ( vtors(3*ntterm+2,ntors), source = 0.0_wp ) 228 do i=1,ntors 229 read (ich,*) tors(1:6,i),vtors(1:2,i), & 230 & (vtors(3*(j-1)+3,i), & 231 & vtors(3*(j-1)+4,i), & 232 & vtors(3*(j-1)+5,i),j=1,tors(5,i)) 233 enddo 234 do i=1,ntors 235 do j=1,tors(5,i) 236 vtors(3*(j-1)+4,i)=vtors(3*(j-1)+4,i)*pi 237 enddo 238 enddo 239 240 if(nhb.gt.0) then 241 allocate( hb(3,nhb), source = 0 ) 242 allocate( vhb(2,nhb), source = 0.0_wp ) 243 read(ich,'(6(3i5,2x))')hb(1:3,1:nhb) 244 do i=1,nhb 245 if(at(hb(3,i)).eq.1)then 246 call hbpara(10.0d0,5.0d0,qff(hb(1,i)),dum1) 247 vhb(1,i)=dum1*scalehb(at(hb(1,i))) 248 call hbpara(10.0d0,5.0d0,qff(hb(2,i)),dum1) 249 vhb(2,i)=dum1*scalehb(at(hb(2,i))) 250 else 251 call hbpara(-6.5d0,1.0d0,qff(hb(3,i)),dum1) 252 vhb(1,i)=scalexb(at(hb(3,i)))*dum1 253 endif 254 enddo 255 endif 256 257 allocate ( nci(3,nnci), source = 0 ) 258 read(ich,'(8(3i5,2x))')nci(1:3,1:nnci) 259 260 close(ich) 261 262 qff=qff*qdum 263 264 write(output_unit,*)'nbond,nangl,ntors :',nbond,nangl,ntors 265 write(output_unit,*)'NCI terms added :',nnci 266 write(output_unit,*)'HB/XB terms added :',nhb 267 return 268 269 do i=1,nbond 270 write(*,'(2i6,3F14.9)')bond(1,i),bond(2,i),vbond(1:3,i) 271 enddo 272 do i=1,nangl 273 write(*,'(3i6,5F14.9)')angl(1,i),angl(2,i),angl(3,i), & 274 & vangl(1:2,i) 275 enddo 276 do i=1,ntors 277 write(*,'(5i6,6F12.8)')tors(1:5,i),vtors(1:2+ntterm,i) 278 enddo 279 280end subroutine rdsolvff 281 282! evaluate the internal energy using the FF as 283! specified in common (ffparam) 284 285subroutine ff_eg(n,at,xyz,e,g) 286 implicit none 287 integer :: n,at(n) 288 real(wp) :: xyz(3,n),e,g(3,n) 289 290! include 'ffparam.f90' 291 292 real(wp) :: vp(3),dt,deddt,rp,cosa,rab2,rcb2,rmul1,rmul2 293 real(wp) :: cosna(4),sinna(4),v(4),sa(4),ca(4),phix(4),dphi(4) 294 real(wp) :: va(3),vb(3),vc(3),vd(3),vba(3),vcb(3),vdc(3),vab(3) 295 real(wp) :: vca(3),vdb(3),vt(3),vu(3),vtu(3),dedt(3),dedu(3) 296 real(wp) :: deda(3),dedb(3),dedc(3),dedd(3),vtmp1(3),vtmp2(3) 297 real(wp) :: rt2,ru2,rtu,rcb,rtru,damp,damp2,rac2,racut,omega 298 real(wp) :: aai,aai2,r2,rjk,dampjk,damp2jk,dampij,damp2ij,dampjl 299 real(wp) :: rkl,rjl,dampkl,damp2kl,damp2jl,term1(3),term2(3),term3(3) 300 integer :: ib,id 301 real(wp) :: dda(3),ddb(3),ddc(3),ddd(3) 302 303 real(wp) :: r,thab,thbc,ra(3),fac,eb,or,kij,rij,alpha,dij 304 real(wp) :: dei(3),dej(3),dek(3),kijk,rb(3),ea,dedtheta,pi6,ban 305 real(wp) :: c0,c1,c2,theta,sinth,costh,expo 306 real(wp) :: dphidri(3),dphidrj(3),dphidrk(3),dphidrl(3),step 307 real(wp) :: et,phi0,rn,kijkl,dedphi,phi,cosrph,cosrph0,er,el 308 real(wp) :: dphi1,dphi2,x1cos,x2cos,x1sin,x2sin,phipi,e1,e2,ef 309 integer :: i,j,k,l,m,ic,ia,list(4),jj,mm,ii,kk,it,nt,it2 310 311 312 e=0.0_wp 313 g=0.0_wp 314 315 ! strech 316 do m=1,nbond 317 i=bond(1,m) 318 j=bond(2,m) 319 r2= ((xyz(1,i)-xyz(1,j))**2 & 320 & +(xyz(2,i)-xyz(2,j))**2 & 321 & +(xyz(3,i)-xyz(3,j))**2) 322 r =sqrt(r2) 323 do ic=1,3 324 ra(ic)=xyz(ic,i)-xyz(ic,j) 325 end do 326 rij=vbond(1,m) 327 kij=vbond(2,m) 328 ! LJ 329 if(bond(3,m).eq.0)then 330 aai=vbond(3,m) 331 aai2 =aai/2 332 e=e+ kij*(1.+(rij/r)**aai - 2.*(rij/r)**aai2 ) 333 fac=aai*kij*( -(rij/r)**aai + (rij/r)**aai2 )/r2 334 ! Morse 335 elseif(bond(3,m).eq.1)then 336 aai=0.5d0*vbond(3,m)/rij 337 e=e+ kij*(1.d0-exp(-aai*(r-rij)))**2 338 fac=2.d0*aai*kij*exp(-aai*(r-rij))*(1.d0-exp(-aai*(r-rij)))/r 339 ! Morse metal 340! elseif(bond(3,m).eq.1)then 341! aai=((7.d0*vbond(3,m)**4+ & 342! & 36.d0*vbond(3,m)**3 & 343! & +44.d0*vbond(3,m)**2) & 344! & /(192.d0*rij**4))**0.25d0 345! e=e+ kij*(1.d0-exp(-aai*(r-rij)))**4 346! fac=4.d0*aai*kij*exp(-aai*(r-rij))* & 347! & (1.d0-exp(-aai*(r-rij)))**3/r 348 endif 349 350 do ic=1,3 351 g(ic,i)=g(ic,i)+fac*ra(ic) 352 g(ic,j)=g(ic,j)-fac*ra(ic) 353 enddo 354 enddo 355 356 ! bend 357 do m=1,nangl 358 j = angl(1,m) 359 i = angl(2,m) 360 k = angl(3,m) 361 c0 =vangl(1,m) 362 kijk=vangl(2,m) 363 va(1:3) = xyz(1:3,i) 364 vb(1:3) = xyz(1:3,j) 365 vc(1:3) = xyz(1:3,k) 366 vab = va-vb 367 vcb = vc-vb 368 rab2 = vab(1)*vab(1) + vab(2)*vab(2) + vab(3)*vab(3) 369 rcb2 = vcb(1)*vcb(1) + vcb(2)*vcb(2) + vcb(3)*vcb(3) 370 call crprod(vcb,vab,vp) 371 rp = norm2(vp)+1.d-14 372 call impsc(vab,vcb,cosa) 373 cosa = dble(min(1.0d0,max(-1.0d0,cosa))) 374 theta= dacos(cosa) 375 376 call abdamp(at(i),at(j),rab2,dampij,damp2ij) 377 call abdamp(at(k),at(j),rcb2,dampjk,damp2jk) 378 damp=dampij*dampjk 379 380 if(pi-c0.lt.1.d-6)then 381 dt = theta - c0 382 ea = kijk * dt**2 383 deddt = 2.d0 * kijk * dt 384 else 385 ea=kijk*(cosa-cos(c0))**2 386 deddt=2.*kijk*sin(theta)*(cos(c0)-cosa) 387 endif 388 389 e = e + ea * damp 390 call crprod(vab,vp,deda) 391 rmul1 = -deddt / (rab2*rp) 392 deda = rmul1*deda 393 call crprod(vcb,vp,dedc) 394 rmul2 = deddt / (rcb2*rp) 395 dedc = rmul2*dedc 396 dedb = deda + dedc 397 term1(1:3)=ea*damp2ij*dampjk*vab(1:3) 398 term2(1:3)=ea*damp2jk*dampij*vcb(1:3) 399 g(1:3,i) = g(1:3,i) + deda(1:3)*damp+term1(1:3) 400 g(1:3,j) = g(1:3,j) - dedb(1:3)*damp-term1(1:3)-term2(1:3) 401 g(1:3,k) = g(1:3,k) + dedc(1:3)*damp+term2(1:3) 402 enddo 403 404 ! torsion 405 do m=1,ntors 406 i=tors(1,m) 407 j=tors(2,m) 408 k=tors(3,m) 409 l=tors(4,m) 410 nt=tors(5,m) 411 phi0 =vtors(1,m) 412 413 if(tors(6,m).ne.2)then 414! ---------------------------------------------------- 415 vab(1:3) = xyz(1:3,i)-xyz(1:3,j) 416 vcb(1:3) = xyz(1:3,j)-xyz(1:3,k) 417 vdc(1:3) = xyz(1:3,k)-xyz(1:3,l) 418 rij = vab(1)*vab(1) + vab(2)*vab(2) + vab(3)*vab(3) 419 rjk = vcb(1)*vcb(1) + vcb(2)*vcb(2) + vcb(3)*vcb(3) 420 rkl = vdc(1)*vdc(1) + vdc(2)*vdc(2) + vdc(3)*vdc(3) 421 call abdamp(at(i),at(j),rij,dampij,damp2ij) 422 call abdamp(at(k),at(j),rjk,dampjk,damp2jk) 423 call abdamp(at(k),at(l),rkl,dampkl,damp2kl) 424 damp= dampjk*dampij*dampkl 425 426 phi=valijklff(n,xyz,i,j,k,l) 427 call dphidr(n,xyz,i,j,k,l,phi,dda,ddb,ddc,ddd) 428 429 et=0 430 dij=0 431 mm=3 432 do it=1,nt 433 rn=vtors(mm,m) 434 dphi1=phi-phi0 435 dphi2=phi+phi0-pi2 436 c1=rn*dphi1+vtors(mm+1,m) 437 c2=rn*dphi2+vtors(mm+1,m) 438 x1cos=cos(c1) 439 x2cos=cos(c2) 440 x1sin=sin(c1) 441 x2sin=sin(c2) 442 phipi=phi-pi 443 ef =erf(phipi) 444 e1 =vtors(mm+2,m)*(1.+x1cos) 445 e2 =vtors(mm+2,m)*(1.+x2cos) 446 et =et+0.5*(1.-ef)*e1+(0.5+0.5*ef)*e2 447 expo =exp(-phipi**2)/spi 448 dij =dij-expo*e1- 0.5*(1.-ef)*vtors(mm+2,m)*x1sin*rn+ & 449 & expo*e2-(0.5+0.5*ef)*vtors(mm+2,m)*x2sin*rn 450 mm=mm+3 451 enddo 452 453 et= et*vtors(2,m) 454 dij=dij*vtors(2,m)*damp 455 456 term1(1:3)=et*damp2ij*dampjk*dampkl*vab(1:3) 457 term2(1:3)=et*damp2jk*dampij*dampkl*vcb(1:3) 458 term3(1:3)=et*damp2kl*dampij*dampjk*vdc(1:3) 459 460 g(1:3,i)=g(1:3,i)+dij*dda(1:3)+term1 461 g(1:3,j)=g(1:3,j)+dij*ddb(1:3)-term1+term2 462 g(1:3,k)=g(1:3,k)+dij*ddc(1:3)+term3-term2 463 g(1:3,l)=g(1:3,l)+dij*ddd(1:3)-term3 464 e=e+et*damp 465! ---------------------------------------------------- 466 else 467! ---------------------------------------------------- 468 vab(1:3) = xyz(1:3,j)-xyz(1:3,i) 469 vcb(1:3) = xyz(1:3,j)-xyz(1:3,k) 470 vdc(1:3) = xyz(1:3,j)-xyz(1:3,l) 471 rij = vab(1)*vab(1) + vab(2)*vab(2) + vab(3)*vab(3) 472 rjk = vcb(1)*vcb(1) + vcb(2)*vcb(2) + vcb(3)*vcb(3) 473 rjl = vdc(1)*vdc(1) + vdc(2)*vdc(2) + vdc(3)*vdc(3) 474 call abdamp(at(i),at(j),rij,dampij,damp2ij) 475 call abdamp(at(k),at(j),rjk,dampjk,damp2jk) 476 call abdamp(at(j),at(l),rjl,dampjl,damp2jl) 477 damp= dampjk*dampij*dampjl 478 479 phi=omega(n,xyz,i,j,k,l) 480 call domegadr(n,xyz,i,j,k,l,phi,dda,ddb,ddc,ddd) 481 482 rn=vtors(3,m) 483 if(rn.gt.1.d-6)then 484 dphi1=phi-phi0 485 c1=dphi1+pi 486 x1cos=cos(c1) 487 x1sin=sin(c1) 488 et =(1.+x1cos)*vtors(2,m) 489 dij =-x1sin*vtors(2,m)*damp 490 else 491 et = vtors(2,m)*(cos(phi) -cos(phi0))**2 492 dij=2.*vtors(2,m)* sin(phi)*(cos(phi0)-cos(phi))*damp 493 endif 494 495 term1(1:3)=et*damp2ij*dampjk*dampjl*vab(1:3) 496 term2(1:3)=et*damp2jk*dampij*dampjl*vcb(1:3) 497 term3(1:3)=et*damp2jl*dampij*dampjk*vdc(1:3) 498 499 g(1:3,i)=g(1:3,i)+dij*dda(1:3)-term1 500 g(1:3,j)=g(1:3,j)+dij*ddb(1:3)+term1+term2+term3 501 g(1:3,k)=g(1:3,k)+dij*ddc(1:3)-term2 502 g(1:3,l)=g(1:3,l)+dij*ddd(1:3)-term3 503 e=e+et*damp 504 endif 505! ---------------------------------------------------- 506 507 enddo 508 509end subroutine ff_eg 510 511 512! nonbonded part 513 514subroutine ff_nonb(n,at,xyz,enb,g) 515 implicit none 516 integer n,at(*) 517 real(wp)xyz(3,n),enb,g(3,n) 518 519! include 'ffparam.f90' 520 521 integer i1,i2,iz1,iz2,k,nk 522 real(wp) r2,r,r4,r6,r06,R0,t6,t8,c6t6,c6t8,t27,drij,e 523 real(wp) dx,dy,dz,c6,x,alpha,oner,aa,damp,damp2,e0 524 525 e=0 526 if(nnci.lt.1)return 527 528 do k=1,nnci 529 i1=nci(1,k) 530 i2=nci(2,k) 531 nk=nci(3,k) 532 dx=xyz(1,i1)-xyz(1,i2) 533 dy=xyz(2,i1)-xyz(2,i2) 534 dz=xyz(3,i1)-xyz(3,i2) 535 r2=dx*dx+dy*dy+dz*dz 536 r =sqrt(r2) 537 oner =1.0d0/r 538 539 iz1=at(i1) 540 iz2=at(i2) 541 R0=r094ff(iz1,iz2) 542 c6=c6ff(i2,i1) 543 r4=r2*r2 544 r6=r4*r2 545 r06=R0**6 546 t6=r6+r06 547 t8=r6*r2+r06*R0*R0 548 c6t6=c6/t6 549 c6t8=c6/t8 550 t27=sr42ff(iz1,iz2)*c6t8 551 e0=c6t6+t27 552 e=e-e0*eps2(nk) 553 drij=eps2(nk)*(c6t6*6.0d0*r4/t6+8.0d0*t27*r6/t8) 554 g(1,i1)=g(1,i1)+dx*drij 555 g(2,i1)=g(2,i1)+dy*drij 556 g(3,i1)=g(3,i1)+dz*drij 557 g(1,i2)=g(1,i2)-dx*drij 558 g(2,i2)=g(2,i2)-dy*drij 559 g(3,i2)=g(3,i2)-dz*drij 560 561 e0=qff(i1)*qff(i2)*oner*eps1(nk) 562 e=e+e0 563 drij=e0/r2 564 g(1,i1)=g(1,i1)-dx*drij 565 g(2,i1)=g(2,i1)-dy*drij 566 g(3,i1)=g(3,i1)-dz*drij 567 g(1,i2)=g(1,i2)+dx*drij 568 g(2,i2)=g(2,i2)+dy*drij 569 g(3,i2)=g(3,i2)+dz*drij 570 571 if(r.lt.25)then 572 x =zabff(iz1,iz2) 573 alpha=r0abff(iz1,iz2) 574 t27 =x*dexp(-alpha*r) 575 e0 =t27*oner 576 e =e + e0*eps2(nk) 577 drij=eps2(nk)*t27*(alpha*r+1.0d0)*oner/r2 578 g(1,i1)=g(1,i1)-dx*drij 579 g(2,i1)=g(2,i1)-dy*drij 580 g(3,i1)=g(3,i1)-dz*drij 581 g(1,i2)=g(1,i2)+dx*drij 582 g(2,i2)=g(2,i2)+dy*drij 583 g(3,i2)=g(3,i2)+dz*drij 584 endif 585 586 enddo 587 588 enb = enb + e 589 590end subroutine ff_nonb 591 592!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 593! HB/XB energy and gradient 594!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 595 596subroutine ff_hb(n,at,xyz,eh,g) 597 implicit none 598 integer n,at(*) 599 real(wp)xyz(3,n),g(3,n),eh 600! include 'ffparam.f90' 601 602 integer i1,i2,i,k,j 603 real(wp)c1,c2,r,er,el,step,edum,dum 604 605 if(nhb.lt.1)return 606 607 do k=1,nhb 608 i1 =hb(1,k) 609 i2 =hb(2,k) 610! r =sqrt((xyz(1,i1)-xyz(1,i2))**2 & 611! & +(xyz(2,i1)-xyz(2,i2))**2 & 612! & +(xyz(3,i1)-xyz(3,i2))**2) 613 r = sqrt(sum((xyz(:,i1)-xyz(:,i2))**2)) 614 if(r.gt.25.0d0)cycle 615 i =hb(3,k) 616 if(at(i).eq.1)then 617 c1 =vhb(1,k) 618 c2 =vhb(2,k) 619 call eabhag(n,i1,i2,i,xyz,c1,c2,eh,g) 620 else 621 c1 =vhb(1,k) 622 call eabxag(n,i1,i2,i,xyz,c1,eh,g) 623 endif 624 enddo 625 626end subroutine ff_hb 627 628!cccccccccccccccccccccccccccccccccccccccccccccc 629! damping of bend and torsion for long 630! bond distances to allow proper dissociation 631!cccccccccccccccccccccccccccccccccccccccccccccc 632 633pure subroutine abdamp(ati,atj,r2,damp,ddamp) 634 use xtb_mctc_convert, only : autoaa 635 use xtb_param_atomicrad, only : atomicRad 636 implicit none 637 integer, intent(in) :: ati,atj 638 real(wp),intent(in) :: r2 639 real(wp),intent(out) :: damp,ddamp 640 real(wp) :: rr,rcut 641 642 ! cut-off at about double of R_cov 643 rcut =3.0 * 3.5710642*((atomicRad(ati)+atomicRad(atj))*autoaa)**2 644 rr =(r2/rcut)**2 645 damp = 1.0d0/(1.0d0+rr) 646 ddamp=-2.d0*2*rr/(r2*(1.0d0+rr)**2) 647 648end subroutine abdamp 649 650!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 651 652subroutine eabh0(n,A,B,H,rab,xyz,ca,cb,energy) 653 implicit none 654 integer A,B,H,n 655 real(wp)xyz(3,n),energy,rab 656 657 real(wp)ang,xy,cosabh,d2ij,d2ik,d2jk,term,temp,fa,fb 658 real(wp)aterm,dampm,dampl,xm,ym,zm,rhm,rab2,ca,cb,da 659 real(wp):: longcut=8. 660 real(wp):: alp= 12 661 real(wp):: alp3= 6 662 663 ! AB distance 664 rab2=rab*rab 665 ! long damping 666 dampl=1./(1.+(rab/longcut)**alp) 667 668 ! cos angle A-B and H (or B-H H) is term 669! D2IK = (XYZ(1,A)-XYZ(1,H))**2+ & 670! & (XYZ(2,A)-XYZ(2,H))**2+ & 671! & (XYZ(3,A)-XYZ(3,H))**2 672! D2JK = (XYZ(1,H)-XYZ(1,B))**2+ & 673! & (XYZ(2,H)-XYZ(2,B))**2+ & 674! & (XYZ(3,H)-XYZ(3,B))**2 675 d2ik = sum((xyz(:,A)-xyz(:,H))**2) 676 d2jk = sum((xyz(:,H)-xyz(:,B))**2) 677 if(d2ik.gt.d2jk)then 678 xy = sqrt(rab2*d2jk) 679 term = 0.5d0 * (rab2+d2jk-d2ik) / xy 680 else 681 xy = sqrt(rab2*d2ik) 682 term = 0.5d0 * (rab2+d2ik-d2jk) / xy 683 endif 684 ! angle damping term 685 aterm = (0.5d0*(term+1.0d0))**alp3 686 687 ! donor-acceptor term 688 fa=d2jk/d2ik 689 fb=d2ik/d2jk 690 da=(ca*fb+cb*fa)/(fa+fb) 691 692 energy=-da*dampl*aterm/(rab2*rab) 693 694end subroutine eabh0 695 696!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 697 698subroutine eabh(n,A,B,H,xyz,ca,cb,energy) 699 implicit none 700 integer A,B,H,n 701 real(wp)xyz(3,n),energy 702 703 real(wp)rab,ang,xy,cosabh,d2ij,d2ik,d2jk,term,temp,fa,fb 704 real(wp)aterm,dampm,dampl,xm,ym,zm,rhm,rab2,ca,cb,da 705 real(wp):: longcut=8. 706 real(wp):: alp=12 707 real(wp):: alp3=6 708 709 ! AB distance 710! rab2= ((xyz(1,A)-xyz(1,B))**2 & 711! & +(xyz(2,A)-xyz(2,B))**2 & 712! & +(xyz(3,A)-xyz(3,B))**2) 713 rab2 = sum((xyz(:,A)-xyz(:,B))**2) ! same as above, but easier to read 714 rab=sqrt(rab2) ! norm2 intrinsic would work here 715 ! long damping 716 dampl=1./(1.+(rab/longcut)**alp) 717 718 ! cos angle A-B and H (or B-H H) is term 719! D2IK = (XYZ(1,A)-XYZ(1,H))**2+ & 720! & (XYZ(2,A)-XYZ(2,H))**2+ & 721! & (XYZ(3,A)-XYZ(3,H))**2 722! D2JK = (XYZ(1,H)-XYZ(1,B))**2+ & 723! & (XYZ(2,H)-XYZ(2,B))**2+ & 724! & (XYZ(3,H)-XYZ(3,B))**2 725 d2ik = sum((xyz(:,A)-xyz(:,H))**2) 726 d2jk = sum((xyz(:,H)-xyz(:,B))**2) 727 if(d2ik.gt.d2jk)then 728 xy = sqrt(rab2*d2jk) 729 term = 0.5d0 * (rab2+d2jk-d2ik) / xy 730 else 731 xy = sqrt(rab2*d2ik) 732 term = 0.5d0 * (rab2+d2ik-d2jk) / xy 733 endif 734 ! angle damping term 735 aterm = (0.5d0*(term+1.0d0))**alp3 736 737 ! donor-acceptor term 738 fa=d2jk/d2ik 739 fb=d2ik/d2jk 740 da=(ca*fb+cb*fa)/(fa+fb) 741 742 ! r^3 only slightly better than r^4 (SG, 8/12) 743 energy=-da*dampl*aterm/(rab2*rab) 744 745end subroutine eabh 746 747!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 748! charge depdendent HB/XB scaling routine 749! for HB a= 10, b=5 750! for XB a=-10, b=1 751!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 752 753subroutine hbpara(a,b,q,c) 754 implicit none 755 real(wp)q,c 756 real(wp)a,b 757 758 c=exp(-a*q)/(exp(-a*q)+b) 759 760end subroutine hbpara 761 762!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 763! hydrogen bonding term with analytical gradient 764!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 765 766subroutine eabhag(n,A,B,H,xyz,ca,cb,energy,gdr) 767 implicit none 768 integer A,B,H,n 769 real(wp)xyz(3,n),ca,cb,energy,gdr(3,n) 770 771 real(wp)cosabh,aterm,rdampl,da 772 real(wp)rab,rah,rbh,rab2,rah2,rbh2,rah4,rbh4 773 real(wp)drah(3),drbh(3),drab(3) 774 real(wp)ga(3),gb(3),gh(3),dg(3),dg2(3) 775 real(wp)gi,denom,ratio,apref,aprod 776 real(wp)eabh 777 778 real(wp):: longcut=8.d0 779 real(wp):: alp= 12.d0 780 real(wp):: alp3= 6.d0 781 782 ! compute distances 783 drah(1:3)=xyz(1:3,A)-xyz(1:3,H) 784 drbh(1:3)=xyz(1:3,B)-xyz(1:3,H) 785 drab(1:3)=xyz(1:3,A)-xyz(1:3,B) 786 787 ! A-B distance 788 rab2=sum(drab**2) 789 rab =sqrt(rab2) 790 791 ! A-H distance 792 rah2=sum(drah**2) 793 rah =sqrt(rah2) 794 795 ! B-H distance 796 rbh2=sum(drbh**2) 797 rbh =sqrt(rbh2) 798 799 ! long damping 800 ratio = (rab/longcut)**alp 801 rdampl=1.d0/(1.d0+ratio) 802 rdampl=rdampl/rab2/rab 803 804 ! cos angle A-B and H (or B-H H) is term 805 if(rah2.gt.rbh2)then 806 aprod = 1.d0/rbh/rab 807 cosabh = -sum(drbh*drab)*aprod 808 else 809 aprod = 1.d0/rah/rab 810 cosabh = sum(drah*drab)*aprod 811 endif 812 813 ! angle damping term 814 aterm = 0.5d0*(cosabh+1.d0) 815 apref = aterm**(alp3-1) 816 aterm = aterm*apref 817 apref = alp3*0.5d0*apref 818 819 ! donor-acceptor term 820 rah4 = rah2*rah2 821 rbh4 = rbh2*rbh2 822 denom = 1.d0/(rah4+rbh4) 823 da = (ca*rah4 + cb*rbh4)*denom 824 825 ! r^3 only slightly better than r^4 (SG, 8/12) 826 eabh = -da*rdampl*aterm 827 828 if(eabh.gt.-1.d-8) return 829 830 energy = energy + eabh 831 832 ! gradient 833 ! donor-acceptor part 834 gi = 4.d0*(ca-cb)*rah2*rbh4*denom*denom 835 gi = -gi*rdampl*aterm 836 ga(1:3) = gi*drah(1:3) 837 gi = 4.d0*(cb-ca)*rbh2*rah4*denom*denom 838 gi = -gi*rdampl*aterm 839 gb(1:3) = gi*drbh(1:3) 840 gh(1:3) = -ga(1:3)-gb(1:3) 841 842 ! long-range damping part 843 gi = rdampl*rdampl*rab*(3.d0+(3.d0+alp)*ratio) 844 gi = gi*da*aterm 845 dg(1:3) = gi*drab(1:3) 846 ga(1:3) = ga(1:3) + dg(1:3) 847 gb(1:3) = gb(1:3) - dg(1:3) 848 849 ! angle part 850 if(rah2.gt.rbh2)then 851 dg(1:3) = -aprod*drbh(1:3) - cosabh*drab(1:3)/rab2 852 dg2(1:3) = aprod*drab(1:3) + cosabh*drbh(1:3)/rbh2 853 gi = -da*rdampl*apref 854 dg(1:3) = gi*dg(1:3) 855 dg2(1:3) = gi*dg2(1:3) 856 857 ga(1:3) = ga(1:3) + dg(1:3) 858 gh(1:3) = gh(1:3) + dg2(1:3) 859 gb(1:3) = gb(1:3) - dg(1:3) - dg2(1:3) 860 861 else 862 dg(1:3) = -aprod*drah(1:3) + cosabh*drab(1:3)/rab2 863 dg2(1:3) = -aprod*drab(1:3) + cosabh*drah(1:3)/rah2 864 gi = -da*rdampl*apref 865 dg(1:3) = gi*dg(1:3) 866 dg2(1:3) = gi*dg2(1:3) 867 868 gb(1:3) = gb(1:3) + dg(1:3) 869 gh(1:3) = gh(1:3) + dg2(1:3) 870 ga(1:3) = ga(1:3) - dg(1:3) - dg2(1:3) 871 872 endif 873 874 ! move gradients into place 875 876 gdr(1:3,A) = gdr(1:3,A) + ga(1:3) 877 gdr(1:3,B) = gdr(1:3,B) + gb(1:3) 878 gdr(1:3,H) = gdr(1:3,H) + gh(1:3) 879 880 return 881end subroutine eabhag 882 883!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 884! XB energy routine 885!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 886 887subroutine eabx(n,A,B,H,xyz,ca,energy) 888 implicit none 889 integer A,B,H,n 890 real(wp)xyz(3,n),energy,rab 891 892 real(wp)ang,xy,cosabh,d2ij,d2ik,d2jk,term,temp,fa,fb 893 real(wp)aterm,dampm,dampl,xm,ym,zm,rhm,rab2,ca,cb,da 894 real(wp):: longcut=120. 895 real(wp):: alp= 6 896 real(wp):: alp3=6 897 898 ! AB distance 899! rab2 = (XYZ(1,A)-XYZ(1,B))**2+ & 900! & (XYZ(2,A)-XYZ(2,B))**2+ & 901! & (XYZ(3,A)-XYZ(3,B))**2 902 rab2 = sum((xyz(:,A)-xyz(:,B))**2) ! same as above, but easier to read 903 ! long damping 904 dampl=1./(1.+(rab2/longcut)**alp) 905 906 ! cos angle A-B and H (or B-H H) is term 907! D2IK = (XYZ(1,A)-XYZ(1,H))**2+ & 908! & (XYZ(2,A)-XYZ(2,H))**2+ & 909! & (XYZ(3,A)-XYZ(3,H))**2 910! D2JK = (XYZ(1,H)-XYZ(1,B))**2+ & 911! & (XYZ(2,H)-XYZ(2,B))**2+ & 912! & (XYZ(3,H)-XYZ(3,B))**2 913 d2ik = sum((xyz(:,A)-xyz(:,H))**2) 914 d2jk = sum((xyz(:,H)-xyz(:,B))**2) 915 if(d2ik.gt.d2jk)then 916 xy = sqrt(rab2*d2jk) 917 term = 0.5d0 * (rab2+d2jk-d2ik) / xy 918 else 919 xy = sqrt(rab2*d2ik) 920 term = 0.5d0 * (rab2+d2ik-d2jk) / xy 921 endif 922 ! angle damping term 923 aterm = (0.5d0*(term+1.0d0))**alp3 924 925 energy=-ca*dampl*aterm/d2jk 926 927end subroutine eabx 928 929!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 930! halogen bonding term with numerical gradient 931!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 932 933subroutine eabxag(n,A,B,H,xyz,ca,energy,gdr) 934 implicit none 935 integer A,B,H,n 936 real(wp)xyz(3,n),ca,energy,gdr(3,n) 937 938 real(wp)cosabh,aterm,rdampl 939 real(wp)rab,rah,rbh,rab2,rah2,rbh2,rah4,rbh4 940 real(wp)drah(3),drbh(3),drab(3) 941 real(wp)ga(3),gb(3),gh(3),dg(3),dg2(3) 942 real(wp)gi,denom,ratio,apref,aprod 943 real(wp)eabh 944 945 real(wp):: longcut=120.d0 946 real(wp):: alp= 6.d0 947 real(wp):: alp3=6.d0 948 949 real(wp)step,er,el,dum 950 integer i,j,i1,i2 951 952 i1=A 953 i2=B 954 i =H 955 call eabx(n,i1,i2,i,xyz,ca,er) 956 energy=energy+er 957 958 step=1.d-6 959 dum=1./(2.0d0*step) 960 do j=1,3 961 xyz(j,i1)=xyz(j,i1)+step 962 call eabx(n,i1,i2,i,xyz,ca,er) 963 xyz(j,i1)=xyz(j,i1)-step*2.0d0 964 call eabx(n,i1,i2,i,xyz,ca,el) 965 xyz(j,i1)=xyz(j,i1)+step 966 gdr(j,i1)=gdr(j,i1)+(er-el)*dum 967 enddo 968 do j=1,3 969 xyz(j,i2)=xyz(j,i2)+step 970 call eabx(n,i1,i2,i,xyz,ca,er) 971 xyz(j,i2)=xyz(j,i2)-step*2.0d0 972 call eabx(n,i1,i2,i,xyz,ca,el) 973 xyz(j,i2)=xyz(j,i2)+step 974 gdr(j,i2)=gdr(j,i2)+(er-el)*dum 975 enddo 976 do j=1,3 977 xyz(j,i )=xyz(j,i )+step 978 call eabx(n,i1,i2,i,xyz,ca,er) 979 xyz(j,i )=xyz(j,i )-step*2.0d0 980 call eabx(n,i1,i2,i,xyz,ca,el) 981 xyz(j,i )=xyz(j,i )+step 982 gdr(j,i )=gdr(j,i )+(er-el)*dum 983 enddo 984 985 return 986end subroutine eabxag 987 988 989!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc 990! set all bonds to Morse function if De > thr 991! ie weak bonds are still LJ 992!cccccccccccccccccccccccccccccccccccccccccccccccccccccccc 993 994subroutine setmorse(n,at,thr,echo) 995 implicit none 996 integer n,at(n) 997 logical echo 998! include 'ffparam.f90' 999 real(wp)thr 1000 integer i,k,i1,i2 1001 1002 k=0 1003 do i=1,nbond12 1004 i1=bond(1,i) 1005 i2=bond(2,i) 1006 if(vbond(2,i).gt.thr) then 1007 bond(3,i)=1 1008 k=k+1 1009 else 1010 bond(3,i)=0 1011 endif 1012 enddo 1013 1014 if(echo)then 1015 write(*,'('' LJ->Morse switch De threshold'',f8.3)') thr 1016 write(*,'('' switching '',i4,'' 1,2 stretch potentials'')') k 1017 endif 1018 1019end subroutine setmorse 1020 1021 1022real(wp)function valijklff(natoms,xyz,i,j,k,l) 1023 1024! ..................................................................... 1025 1026implicit none 1027 1028external & 1029 & vecnorm,valijk 1030 1031integer & 1032 & ic,i,j,k,l,natoms 1033 1034real(wp)& 1035 & xyz(3,natoms), & 1036 & eps,ra(3),rb(3),rc(3),na(3),nb(3), & 1037 & rab,rbc,thab,thbc,valijk, & 1038 & vecnorm,nan,nbn,rcn,snanb,deter,pi 1039 1040parameter (eps=1.0d-14) 1041data pi/3.1415926535897932384626433832795029d0/ 1042 1043! ... get torsion coordinate 1044do ic=1,3 1045 ra(ic)=xyz(ic,j)-xyz(ic,i) 1046 rb(ic)=xyz(ic,k)-xyz(ic,j) 1047 rc(ic)=xyz(ic,l)-xyz(ic,k) 1048end do 1049 1050! ... determinante of rb,ra,rc 1051deter= ra(1)*(rb(2)*rc(3)-rb(3)*rc(2)) & 1052 & -ra(2)*(rb(1)*rc(3)-rb(3)*rc(1)) & 1053 & +ra(3)*(rb(1)*rc(2)-rb(2)*rc(1)) 1054 1055thab=valijk(natoms,xyz,i,k,j) 1056thbc=valijk(natoms,xyz,j,l,k) 1057call crossprod(ra,rb,na) 1058call crossprod(rb,rc,nb) 1059nan=vecnorm(na,3,1) 1060nbn=vecnorm(nb,3,1) 1061 1062snanb=0.0d0 1063do ic=1,3 1064 snanb=snanb+na(ic)*nb(ic) 1065end do 1066if (abs(abs(snanb)-1.d0).lt.eps) then 1067 snanb=sign(1.d0,snanb) 1068end if 1069 1070valijklff=acos(snanb) 1071 1072! the gradient dphir is only compatible with this subroutine 1073! if the statement below is commented out. If not, opt. and 1074! Hessian show large errors and imags. I don't understand 1075! this entirely but thats how it is. 1076! SG, Sat May 24 11:41:42 CEST 2014 1077 1078! if (deter.lt.0) then 1079! valijkl=2.d0*pi-valijkl 1080! end if 1081 1082end function valijklff 1083 1084subroutine valelff(z) 1085 real(wp)at,z(*) 1086 integer i 1087 1088 z(1:max_elem)=0 1089 do i=1,54 1090 at=float(i) 1091 ! -----H 1092 z(i)=at 1093 if(i.le.10.and.i.gt.2)then 1094 ! ------Li-F 1095 z(i)=at-2 1096 endif 1097 if(i.le.13.and.i.gt.10)then 1098 ! ------Na-Al 1099 z(i)=at-10 1100 endif 1101 if(i.le.18.and.i.gt.13)then 1102 ! ------Si-Ar 1103 z(i)=at-10 1104 endif 1105 if((i.le.36.and.i.ge.30).or.i.eq.19.or. i.eq.20) then 1106 ! ------K,Ca,Zn-Kr 1107 ! 4s4p 1108 z(i)=at-18 1109 if(i.ge.30)z(i)=at-28 1110 endif 1111 if(i.le.29.and.i.ge.21) then 1112 ! Sc-Cu 1113 z(i)=at-18 1114 endif 1115 if(i.le.47.and.i.ge.37) then 1116 ! Y-Ag 1117 z(i)=at-36 1118 endif 1119 if(i.gt.47.and.i.le.54) then 1120 ! In-Xe 1121 z(i)=at-46 1122 endif 1123 enddo 1124 1125 ! not well define for these metals 1126 do i=57,80 1127 z(i)=4.0 1128 enddo 1129 1130 do i=81,86 1131 z(i)=float(i)-78.0 1132 enddo 1133 1134 ! modifications (fit) 1135 z(1:2) =z(1:2)* 2.35 1136 z(3:10) =z(3:10)* 0.95 1137 z(11:18)=z(11:18)*0.75 1138 z(19:54)=z(19:54)*0.65 1139 ! just extrapolated 1140 z(55:max_elem)=z(55:max_elem)*0.60 1141 1142 ! special for group 1 and 2, fitted to E_int and O-M distances 1143 ! at TPSS-D3/def2-TZVP level 1144 z(3) =1.7 1145 z(11)=2.5 1146 z(19)=3.0 1147 z(37)=3.0 1148 z(55)=3.0 1149 z(87)=3.0 1150 1151 z( 4)=5.5 1152 z(12)=3.0 1153 z(20)=2.8 1154 z(38)=2.6 1155 z(56)=2.4 1156 z(88)=2.4 1157 1158end subroutine valelff 1159end module xtb_qmdff 1160