1! This file is part of xtb. 2! 3! Copyright (C) 2019-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/>. 17module xtb_gfnff_ini2 18 use xtb_gfnff_data, only : TGFFData 19 use xtb_gfnff_neighbourlist, only : TGFFNeighbourList 20 use xtb_gfnff_topology, only : TGFFTopology 21 use xtb_type_environment, only : TEnvironment 22 implicit none 23 private 24 public :: gfnff_neigh, getnb, nbondmat 25 public :: pairsbond, pilist, nofs, xatom, ctype, amide, amideH, alphaCO 26 public :: ringsatom, ringsbond, ringsbend, ringstors, ringstorl 27 public :: chktors, chkrng, hbonds, getring36, ssort, goedeckera, qheavy 28 public :: gfnff_hbset, gfnff_hbset0, bond_hbset, bond_hbset0 29 public :: bond_hb_AHB_set, bond_hb_AHB_set1, bond_hb_AHB_set0 30 31contains 32 33subroutine gfnff_neigh(env,makeneighbor,natoms,at,xyz,rab,fq,f_in,f2_in,lintr,mchar,hyb,itag,nbm,nbf,param,topo) 34 use xtb_gfnff_param 35 implicit none 36 character(len=*), parameter :: source = 'gfnff_ini2_neigh' 37 type(TEnvironment), intent(inout) :: env 38 type(TGFFData), intent(in) :: param 39 type(TGFFTopology), intent(inout) :: topo 40 logical makeneighbor 41 integer at(natoms),natoms 42 integer hyb (natoms) 43 integer itag(natoms) 44 integer nbm(20,natoms) ! needed for ring assignment (done without metals) 45 integer nbf(20,natoms) ! full needed for fragment assignment 46 real*8 rab (natoms*(natoms+1)/2) 47 real*8 xyz(3,natoms) 48 real*8 mchar(natoms) 49 real*8 fq 50 real*8 f_in,f2_in ! radius scaling for atoms/metal atoms recpectively 51 real*8 lintr ! threshold for linearity 52 53 logical etacoord,da,strange_iat,metal_iat 54 integer,allocatable :: nbdum(:,:) 55 real*8 ,allocatable :: cn(:),rtmp(:) 56 integer iat,i,j,k,ni,ii,jj,kk,ll,lin,ati,nb20i,nbdiff,hc_crit,nbmdiff,nnf,nni,nh,nm 57 integer ai,aj,nn,im,ncm,l,no 58 real*8 r,pi,a1,f,f1,phi,f2,rco,fat(86) 59 data pi/3.1415926535897932384626433832795029d0/ 60 data fat / 86 * 1.0d0 / 61 62! special hacks 63 fat( 1)=1.02 64 fat( 4)=1.03 65 fat( 5)=1.02 66 fat( 8)=1.02 67 fat( 9)=1.05 68 fat(10)=1.10 69 fat(11)=1.01 70 fat(12)=1.02 71 fat(15)=0.97 72 fat(18)=1.10 73 fat(19)=1.02 74 fat(20)=1.02 75 fat(38)=1.02 76 fat(34)=0.99 77 fat(50)=1.01 78 fat(51)=0.99 79 fat(52)=0.95 80 fat(53)=0.98 81 fat(56)=1.02 82 fat(76)=1.02 83 fat(82)=1.06 84 fat(83)=0.95 85 86 allocate(cn(natoms),rtmp(natoms*(natoms+1)/2),nbdum(20,natoms)) 87 88! determine the neighbor list 89 if(makeneighbor) then 90 91 topo%nb =0 ! without highly coordinates atoms 92 nbm=0 ! without any metal 93 nbf=0 ! full 94 95 do i=1,natoms 96 cn(i)=dble(param%normcn(at(i))) 97 enddo 98 call gfnffrab(natoms,at,cn,rtmp) ! guess RAB based on "normal" CN 99 do i=1,natoms 100 ai=at(i) 101 f1=fq 102 if(param%metal(ai) > 0) f1 = f1 * 2.0d0 103 do j=1,i-1 104 f2=fq 105 aj=at(j) 106 if(param%metal(aj) > 0) f2 = f2 * 2.0d0 107 k=lin(j,i) 108 rco=rtmp(k) 109 rtmp(k)=rtmp(k)-topo%qa(i)*f1-topo%qa(j)*f2 ! change radius of atom i and j with charge 110! element specials 111 rtmp(k)=rtmp(k)*fat(ai)*fat(aj) 112 enddo 113 enddo 114 115 call getnb(natoms,at,rtmp,rab,mchar,1,f_in,f2_in,nbdum,nbf,param) ! full 116 call getnb(natoms,at,rtmp,rab,mchar,2,f_in,f2_in,nbf ,topo%nb,param) ! no highly coordinates atoms 117 call getnb(natoms,at,rtmp,rab,mchar,3,f_in,f2_in,nbf ,nbm,param) ! no metals and unusually coordinated stuff 118 119! take the input 120 else 121 122 nbf = topo%nb 123 nbm = topo%nb 124 125 endif 126! done 127 128 itag = 0 ! save special hyb info 129 130! tag atoms in nb(19,i) if they belong to a cluster (which avoids the ring search) 131 do i=1,natoms 132 if(nbf(20,i).eq.0.and.param%group(at(i)).ne.8)then 133 write(env%unit,'(''!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'')') 134 write(env%unit,'('' warning: no bond partners for atom'',i4)')i 135 write(env%unit,'(''!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'')') 136 endif 137 if(at(i).lt.11.and.nbf(20,i).gt.2)then 138 do k=1,nbf(20,i) 139 kk=nbf(k,i) 140 if(param%metal(at(kk)).ne.0.or.topo%nb(20,kk).gt.4) then 141 topo%nb (19,i)=1 142 nbf(19,i)=1 143 nbm(19,i)=1 144 endif 145 enddo 146 endif 147! write(env%unit,*) i,(topo%nb(j,i),j=1,topo%nb(20,i)) 148 enddo 149 150!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 151! hybridization states 152!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 153 do i=1,natoms 154 ati = at(i) 155! important: determine cases where atom is pi bonded to a metal and thus 156! the hyb must be obtained from the reduced (wo metals) neighbor list 157 etacoord=.false. 158 if(ati.le.10)then 159 if(ati.eq.6.and.nbf(20,i).ge.4.and.nbm(20,i).eq.3) etacoord=.true. ! CP case 160 if(ati.eq.6.and.nbf(20,i).eq.3.and.nbm(20,i).eq.2) etacoord=.true. ! alkyne case 161 nm=0 162 do k=1,nbf(20,i) ! how many metals ? and which 163 kk=nbf(k,i) 164 if(param%metal(at(kk)).ne.0) then 165 nm=nm+1 166 im=kk 167 endif 168 enddo 169 if(nm.eq.0) then 170 etacoord=.false. ! etacoord makes no sense without metals! 171 elseif(nm.eq.1)then ! distinguish M-CR2-R i.e. not an eta coord. 172 ncm=0 173 do k=1,nbf(20,i) ! 174 if(nbf(k,i).ne.im)then ! all neighbors that are not the metal im 175 kk=nbf(k,i) 176 do l=1,nbf(20,kk) 177 if(nbf(l,kk).eq.im) ncm=ncm+1 ! ncm=1 is alkyne, =2 is cp 178 enddo 179 endif 180 enddo 181 if(ncm.eq.0) etacoord=.false. 182 endif 183 endif 184 if(etacoord)then 185 itag(i)=-1 186 nbdum(1:20,i)=nbm(1:20,i) 187 else 188 nbdum(1:20,i)=nbf(1:20,i) ! take full set of neighbors by default 189 endif 190 enddo 191 192 do i=1,natoms 193 ati = at(i) 194 hyb(i)=0 ! don't know it 195 nbdiff =nbf(20,i)-topo%nb (20,i) 196 nbmdiff=nbf(20,i)-nbm(20,i) 197 nb20i=nbdum(20,i) 198 nh=0 199 no=0 200 do j=1,nb20i 201 if(at(nbdum(j,i)).eq.1) nh=nh+1 202 if(at(nbdum(j,i)).eq.8) no=no+1 203 enddo 204! H 205 if(param%group(ati).eq.1) then 206 if(nb20i.eq.2) hyb(i)=1 ! bridging H 207 if(nb20i.gt.2) hyb(i)=3 ! M+ tetra coord 208 if(nb20i.gt.4) hyb(i)=0 ! M+ HC 209 endif 210! Be 211 if(param%group(ati).eq.2) then 212 if(nb20i.eq.2) hyb(i)=1 ! bridging M 213 if(nb20i.gt.2) hyb(i)=3 ! M+ tetra coord 214 if(nb20i.gt.4) hyb(i)=0 ! 215 endif 216! B 217 if(param%group(ati).eq.3) then 218 if(nb20i.gt.4) hyb(i)=3 219 if(nb20i.gt.4.and.ati.gt.10.and.nbdiff.eq.0) hyb(i)=5 220 if(nb20i.eq.4) hyb(i)=3 221 if(nb20i.eq.3) hyb(i)=2 222 if(nb20i.eq.2) hyb(i)=1 223 endif 224! C 225 if(param%group(ati).eq.4) then 226 if(nb20i.ge.4) hyb(i)=3 227 if(nb20i.gt.4.and.ati.gt.10.and.nbdiff.eq.0) hyb(i)=5 228 if(nb20i.eq.3) hyb(i)=2 229 if(nb20i.eq.2) then 230 call bangl(xyz,nbdum(1,i),i,nbdum(2,i),phi) 231 if(phi*180./pi.lt.150.0)then ! geometry dep. setup! GEODEP 232 hyb(i)=2 ! otherwise, carbenes will not be recognized 233 itag(i)=1 ! tag for Hueckel and HB routines 234 else 235 hyb(i)=1 ! linear triple bond etc 236 endif 237 if(topo%qa(i).lt.-0.4) then 238 hyb(i)=2 239 itag(i)=0 ! tag for Hueckel and HB routines 240 endif 241 endif 242 if(nb20i.eq.1) hyb(i)=1 ! CO 243 endif 244! N 245 if(param%group(ati).eq.5) then 246 if(nb20i.ge.4) hyb(i)=3 247 if(nb20i.gt.4.and.ati.gt.10.and.nbdiff.eq.0) hyb(i)=5 248 if(nb20i.eq.3) hyb(i)=3 249 if(nb20i.eq.3.and.ati.eq.7) then 250 kk=0 251 ll=0 252 nn=0 253 do j=1,3 254 jj=nbdum(j,i) 255 if(at(jj).eq. 8.and.topo%nb(20,jj).eq.1) kk=kk+1 ! check for NO2 or R2-N=O 256 if(at(jj).eq. 5.and.topo%nb(20,jj).eq.4) ll=ll+1 ! check for B-N, if the CN(B)=4 the N is loosely bound and sp2 257 if(at(jj).eq.16.and.topo%nb(20,jj).eq.4) nn=nn+1 ! check for N-SO2- 258 enddo 259 if(nn.eq.1.and.ll.eq.0.and.kk.eq.0) hyb(i)=3 260 if(ll.eq.1.and.nn.eq.0) hyb(i)=2 261 if(kk.ge.1) then 262 hyb(i)=2 263 itag(i)=1 ! tag for Hueckel with no el. for the N in NO2 264 endif 265 if(nbmdiff.gt.0.and.nn.eq.0) hyb(i)=2 ! pyridin N coord. to heavy atom 266 endif 267 if(nb20i.eq.2) then 268 hyb(i)=2 269 call bangl(xyz,nbdum(1,i),i,nbdum(2,i),phi) 270 jj=nbdum(1,i) 271 kk=nbdum(2,i) 272 if(nbdum(20,jj).eq.1.and.at(jj).eq.6) hyb(i)=1 ! R-N=C 273 if(nbdum(20,kk).eq.1.and.at(kk).eq.6) hyb(i)=1 ! R-N=C 274 if(nbdum(20,jj).eq.1.and.at(jj).eq.7) hyb(i)=1 ! R-N=N in e.g. diazomethane 275 if(nbdum(20,kk).eq.1.and.at(kk).eq.7) hyb(i)=1 ! R-N=N in e.g. diazomethane 276 if(nbdum(1,i).gt.0.and.param%metal(at(nbdum(1,i))).gt.0) hyb(i)=1 ! M-NC-R in e.g. nitriles 277 if(nbdum(2,i).gt.0.and.param%metal(at(nbdum(2,i))).gt.0) hyb(i)=1 ! M-NC-R in e.g. nitriles 278 if(at(jj).eq.7.and.at(kk).eq.7.and. & 279 & nbdum(20,jj).le.2.and.nbdum(20,kk).le.2) hyb(i)=1 ! N=N=N 280 if(phi*180./pi.gt.lintr) hyb(i)=1 ! geometry dep. setup! GEODEP 281 endif 282 if(nb20i.eq.1) hyb(i)=1 283 endif 284! O 285 if(param%group(ati).eq.6) then 286 if(nb20i.ge.3) hyb(i)=3 287 if(nb20i.gt.3.and.ati.gt.10.and.nbdiff.eq.0) hyb(i)=5 288 if(nb20i.eq.2) hyb(i)=3 289 if(nb20i.eq.2.and.nbmdiff.gt.0) then 290 call nn_nearest_noM(i,natoms,at,topo%nb,rab,j,param) ! CN of closest non-M atom 291 if(j.eq.3) hyb(i)=2 ! M-O-X konj 292 if(j.eq.4) hyb(i)=3 ! M-O-X non 293 endif 294 if(nb20i.eq.1) hyb(i)=2 295 if(nb20i.eq.1.and.nbdiff.eq.0) then 296 if(topo%nb(20,topo%nb(1,i)).eq.1) hyb(i)=1 ! CO 297 endif 298 endif 299! F 300 if(param%group(ati).eq.7) then 301 if(nb20i.eq.2) hyb(i)=1 302 if(nb20i.gt.2.and.ati.gt.10) hyb(i)=5 303 endif 304! Ne 305 if(param%group(ati).eq.8) then 306 hyb(i)=0 307 if(nb20i.gt.0.and.ati.gt.2) hyb(i)=5 308 endif 309! done with main groups 310 if(param%group(ati).le.0) then ! TMs 311 nni=nb20i 312 if(nh.ne.0.and.nh.ne.nni) nni=nni-nh ! don't count Hs 313 if(nni.le.2) hyb(i)=1 314 if(nni.le.2.and.param%group(ati).le.-6) hyb(i)=2 315 if(nni.eq.3) hyb(i)=2 316 if(nni.eq.4.and.param%group(ati).gt.-7) hyb(i)=3 ! early TM, tetrahedral 317 if(nni.eq.4.and.param%group(ati).le.-7) hyb(i)=3 ! late TM, square planar 318 if(nni.eq.5.and.param%group(ati).eq.-3) hyb(i)=3 ! Sc-La are tetrahedral CN=5 319 endif 320 enddo 321 322 topo%nb = nbdum ! list is complete but hyb determination is based only on reduced (without metals) list 323 324 deallocate(nbdum) 325 326 j = 0 327 do i=1,natoms 328 if(topo%nb(20,i).gt.12) j = j +1 329 do k=1,topo%nb(20,i) 330 kk=topo%nb(k,i) 331 if(at(kk).eq.6.and.at(i).eq.6.and.itag(i).eq.1.and.itag(kk).eq.1) then ! check the very special situation of 332 itag(i) =0 ! two carbene C bonded which is an arine 333 itag(kk)=0 334 endif 335 enddo 336 enddo 337 if(dble(j)/dble(natoms).gt.0.3) then 338 call env%error(' too many atoms with extreme high CN', source) 339 end if 340 341 end subroutine gfnff_neigh 342 343!ccccccccccccccccccccccccccccccccccccccccccccccccc 344! fill neighbor list 345!ccccccccccccccccccccccccccccccccccccccccccccccccc 346 347 subroutine getnb(n,at,rad,r,mchar,icase,f,f2,nbf,nb,param) 348 implicit none 349 type(TGFFData), intent(in) :: param 350 integer n,at(n),nbf(20,n),nb(20,n) 351 real*8 rad(n*(n+1)/2),r(n*(n+1)/2),mchar(n),f,f2 352 353 integer i,j,k,nn,icase,hc_crit,nnfi,nnfj,lin 354 integer tag(n*(n+1)/2) 355 real*8 rco,fm 356 357 nb = 0 ! resulting array (nbf is full from first call) 358 tag = 0 359 do i=1,n 360 nnfi=nbf(20,i) ! full CN of i, only valid for icase > 1 361 do j=1,i-1 362 nnfj=nbf(20,j) ! full CN of i 363 fm=1.0d0 364! full case 365 if(icase.eq.1)then 366 if(param%metal(at(i)).eq.2) fm=fm*f2 !change radius for metal atoms 367 if(param%metal(at(j)).eq.2) fm=fm*f2 368 if(param%metal(at(i)).eq.1) fm=fm*(f2+0.025) 369 if(param%metal(at(j)).eq.1) fm=fm*(f2+0.025) 370 endif 371! no HC atoms 372 if(icase.eq.2)then 373 hc_crit = 6 374 if(param%group(at(i)).le.2) hc_crit = 4 375 if(nnfi.gt.hc_crit) cycle 376 hc_crit = 6 377 if(param%group(at(j)).le.2) hc_crit = 4 378 if(nnfj.gt.hc_crit) cycle 379 endif 380! no metals and unusually coordinated stuff 381 if(icase.eq.3)then 382 if(mchar(i).gt.0.25 .or. param%metal(at(i)).gt.0) cycle ! metal case TMonly ?? TODO 383 if(mchar(j).gt.0.25 .or. param%metal(at(j)).gt.0) cycle ! metal case 384 if(nnfi.gt.param%normcn(at(i)).and.at(i).gt.10) cycle ! HC case 385 if(nnfj.gt.param%normcn(at(j)).and.at(j).gt.10) cycle ! HC case 386 endif 387 k=lin(j,i) 388 rco=rad(k) !(rad(i)+rad(j))/0.5291670d0 389! R est. R0 390 if(r(k).lt. fm * f * rco) tag(k)=1 391 enddo 392 enddo 393 394 do i=1,n 395 nn = 0 396 do j=1,n 397 if(tag(lin(j,i)).eq.1.and.i.ne.j) then 398 nn=nn+1 399 nb(nn,i)=j 400 endif 401 enddo 402 nb(20,i) = nn 403 enddo 404 405 end subroutine getnb 406 407!ccccccccccccccccccccccccccccccccccccccccccccccccc 408! find the CN of nearest non metal of atom i 409!ccccccccccccccccccccccccccccccccccccccccccccccccc 410 411 subroutine nn_nearest_noM(ii,n,at,nb,r,nn,param) 412 implicit none 413 type(TGFFData), intent(in) :: param 414 integer ii,n,at(n),nn,nb(20,n) 415 real*8 r(n*(n+1)/2) 416 417 integer jmin,j,jj,lin 418 real*8 rmin 419 420 nn=0 421 rmin=1.d+42 422 jmin=0 423 do j=1,nb(20,ii) 424 jj=nb(j,ii) 425 if(param%metal(at(jj)).ne.0) cycle 426 if(r(lin(jj,ii)).lt.rmin)then 427 rmin=r(lin(jj,ii)) 428 jmin=jj 429 endif 430 enddo 431 432 if(jmin.gt.0) nn=nb(20,jmin) 433 434 end subroutine nn_nearest_noM 435 436!ccccccccccccccccccccccccccccccccccccccccccccccccc 437! find smallest ring in which atom i is located 438!ccccccccccccccccccccccccccccccccccccccccccccccccc 439 440 subroutine ringsatom(n,i,c,s,rings) 441 implicit none 442 integer n,i,k,l,c(10,20,n),s(20,n),rings,rings1 443 444 rings=99 445 do k=1,s(20,i) ! all rings of atom i 446 if(s(k,i).lt.rings) rings=s(k,i) 447 enddo 448 449 end subroutine ringsatom 450 451!ccccccccccccccccccccccccccccccccccccccccccccccccc 452! find smallest ring in which bond i-j is located 453!ccccccccccccccccccccccccccccccccccccccccccccccccc 454 455 subroutine ringsbond(n,i,j,c,s,rings) 456 implicit none 457 integer n,i,j,k,l,c(10,20,n),s(20,n),rings,rings1,rings2 458 459 rings1=99 460 rings2=99 461 do k=1,s(20,i) ! all rings of atom i 462 do l=1,s(k,i) ! all atoms of ring k 463 if(c(l,k,i).eq.j.and.s(k,i).lt.rings1)then 464 rings1=s(k,i) 465 endif 466 enddo 467 enddo 468 do k=1,s(20,j) ! all rings of atom i 469 do l=1,s(k,j) ! all atoms of ring k 470 if(c(l,k,j).eq.i.and.s(k,j).lt.rings2)then 471 rings2=s(k,j) 472 endif 473 enddo 474 enddo 475 continue 476 rings=min(rings1,rings2) 477 if(rings.eq.99) rings=0 478 479 end subroutine ringsbond 480 481!cccccccccccccccccccccccccccccccccccccccccccccccccccc 482! find smallest ring in which angle i-j-k is located 483!cccccccccccccccccccccccccccccccccccccccccccccccccccc 484 485 subroutine ringsbend(n,i,j,k,c,s,rings) 486 implicit none 487 integer n,i,j,k,rings 488 integer c(10,20,n),s(20,n) 489 integer itest,rings1,rings2,rings3,m,l 490 491 if(s(20,i).eq.0.or.s(20,j).eq.0.or.s(20,k).eq.0)then 492 rings=0 493 return 494 endif 495 496 rings1=99 497 rings2=99 498 rings3=99 499 500 do m=1,s(20,i) ! all rings of atom i 501 itest=0 502 do l=1,s(m,i) ! all atoms of ring m 503 if(c(l,m,i).eq.j.or.c(l,m,i).eq.k) itest=itest+1 504 enddo 505 if(itest.eq.2.and.s(m,i).lt.rings1) rings1=s(m,i) 506 enddo 507 do m=1,s(20,j) ! all rings of atom j 508 itest=0 509 do l=1,s(m,j) ! all atoms of ring m 510 if(c(l,m,j).eq.i.or.c(l,m,j).eq.k) itest=itest+1 511 enddo 512 if(itest.eq.2.and.s(m,j).lt.rings2) rings2=s(m,j) 513 enddo 514 do m=1,s(20,k) ! all rings of atom k 515 itest=0 516 do l=1,s(m,k) ! all atoms of ring m 517 if(c(l,m,k).eq.i.or.c(l,m,k).eq.j) itest=itest+1 518 enddo 519 if(itest.eq.2.and.s(m,j).lt.rings3) rings3=s(m,k) 520 enddo 521 522 rings=min(rings1,rings2,rings3) 523 if(rings.eq.99) rings=0 524 525 end subroutine ringsbend 526 527!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 528! find smallest torsion in which angle i-j-k-l is located 529!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 530 531 subroutine ringstors(n,i,j,k,l,c,s,rings) 532 implicit none 533 integer n,i,j,k,l,rings 534 integer c(10,20,n),s(20,n) 535 integer itest,rings1,rings2,rings3,rings4,m,a 536! integer ht1,ht2,ht3,ht4 537 538 if(s(20,i).eq.0.or.s(20,j).eq.0.or.s(20,k).eq.0.or.s(20,l).eq.0)then 539 rings=0 540 return 541 endif 542 543 rings1=99 544 rings2=99 545 rings3=99 546 rings4=99 547 548 do m=1,s(20,i) ! all rings of atom i 549 itest=0 550 do a=1,s(m,i) ! all atoms of ring m 551 if(c(a,m,i).eq.j.or.c(a,m,i).eq.k.or.c(a,m,i).eq.l) itest=itest+1 552 enddo 553 if(itest.eq.3.and.s(m,i).lt.rings1) then 554 rings1=s(m,i) 555! ht1=c(m,19,i) 556 endif 557 enddo 558 do m=1,s(20,j) ! all rings of atom j 559 itest=0 560 do a=1,s(m,j) ! all atoms of ring m 561 if(c(a,m,j).eq.i.or.c(a,m,j).eq.k.or.c(a,m,j).eq.l) itest=itest+1 562 enddo 563 if(itest.eq.3.and.s(m,j).lt.rings2) then 564 rings2=s(m,j) 565! ht2=c(m,19,j) 566 endif 567 enddo 568 do m=1,s(20,k) ! all rings of atom k 569 itest=0 570 do a=1,s(m,k) ! all atoms of ring m 571 if(c(a,m,k).eq.i.or.c(a,m,k).eq.j.or.c(a,m,k).eq.l) itest=itest+1 572 enddo 573 if(itest.eq.3.and.s(m,k).lt.rings3) then 574 rings3=s(m,k) 575! ht3=c(m,19,k) 576 endif 577 enddo 578 do m=1,s(20,l) ! all rings of atom k 579 itest=0 580 do a=1,s(m,l) ! all atoms of ring m 581 if(c(a,m,l).eq.i.or.c(a,m,l).eq.k.or.c(a,m,l).eq.j) itest=itest+1 582 enddo 583 if(itest.eq.3.and.s(m,l).lt.rings4) then 584 rings4=s(m,l) 585! ht4=c(m,19,l) 586 endif 587 enddo 588 589 rings=min(rings1,rings2,rings3,rings4) 590 if(rings.eq.99) then 591 rings=0 592 endif 593 594 end subroutine ringstors 595 596!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 597! find smallest torsion in which angle i-j-k-l is located 598!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 599 600 subroutine ringstorl(n,i,j,k,l,c,s,ringl) 601 implicit none 602 integer n,i,j,k,l,ringl 603 integer c(10,20,n),s(20,n) 604 integer itest,rings1,rings2,rings3,rings4,m,a 605 606 if(s(20,i).eq.0.or.s(20,j).eq.0.or.s(20,k).eq.0.or.s(20,l).eq.0)then 607 ringl=0 608 return 609 endif 610 611 rings1=-99 612 rings2=-99 613 rings3=-99 614 rings4=-99 615 616 do m=1,s(20,i) ! all rings of atom i 617 itest=0 618 do a=1,s(m,i) ! all atoms of ring m 619 if(c(a,m,i).eq.j.or.c(a,m,i).eq.k.or.c(a,m,i).eq.l) itest=itest+1 620 enddo 621 if(itest.eq.3.and.s(m,i).gt.rings1) then 622 rings1=s(m,i) 623 endif 624 enddo 625 do m=1,s(20,j) ! all rings of atom j 626 itest=0 627 do a=1,s(m,j) ! all atoms of ring m 628 if(c(a,m,j).eq.i.or.c(a,m,j).eq.k.or.c(a,m,j).eq.l) itest=itest+1 629 enddo 630 if(itest.eq.3.and.s(m,j).gt.rings2) then 631 rings2=s(m,j) 632 endif 633 enddo 634 do m=1,s(20,k) ! all rings of atom k 635 itest=0 636 do a=1,s(m,k) ! all atoms of ring m 637 if(c(a,m,k).eq.i.or.c(a,m,k).eq.j.or.c(a,m,k).eq.l) itest=itest+1 638 enddo 639 if(itest.eq.3.and.s(m,k).gt.rings3) then 640 rings3=s(m,k) 641 endif 642 enddo 643 do m=1,s(20,l) ! all rings of atom k 644 itest=0 645 do a=1,s(m,l) ! all atoms of ring m 646 if(c(a,m,l).eq.i.or.c(a,m,l).eq.k.or.c(a,m,l).eq.j) itest=itest+1 647 enddo 648 if(itest.eq.3.and.s(m,l).gt.rings4) then 649 rings4=s(m,l) 650 endif 651 enddo 652 653 ringl=max(rings1,rings2,rings3,rings4) 654 if(ringl.eq.-99) then 655 ringl=0 656 endif 657 658 end subroutine ringstorl 659 660!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 661 662 logical function chktors(n,xyz,i,j,k,l) ! true if dihedral angle is bad i.e. near 180 663 implicit none 664 integer n,i,j,k,l 665 real*8 xyz(3,n),phi 666 667 chktors=.true. 668 669 call bangl(xyz,j,i,k,phi) 670! write(env%unit,*) phi*180./3.1415926d0 671 if(phi*180./3.1415926d0.gt.170.0d0) return 672 call bangl(xyz,i,j,l,phi) 673! write(env%unit,*) phi*180./3.1415926d0 674 if(phi*180./3.1415926d0.gt.170.0d0) return 675 676 chktors=.false. 677 678 end function chktors 679 680 logical function chkrng(nn,n,c) 681 implicit none 682 integer n,idum(nn),nn,c(10),i,j 683 chkrng=.true. 684 idum=0 685 do i=1,n 686 idum(c(i))=idum(c(i))+1 687 enddo 688 j=0 689 do i=1,nn 690 if(idum(i).eq.1) j=j+1 691 enddo 692 if(j.ne.n) chkrng=.false. 693 end function chkrng 694 695!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 696 697subroutine gfnff_hbset(n,at,xyz,sqrab,topo,nlist,hbthr1,hbthr2) 698use xtb_mctc_accuracy, only : wp 699 use xtb_gfnff_param 700 implicit none 701 type(TGFFTopology), intent(in) :: topo 702 type(TGFFNeighbourList), intent(inout) :: nlist 703 integer n 704 integer at(n) 705 real(wp) sqrab(n*(n+1)/2) 706 real(wp) xyz(3,n) 707 real(wp), intent(in) :: hbthr1, hbthr2 708 709 integer i,j,k,nh,ia,ix,lin,ij,inh,jnh 710 real(wp) rab,rmsd 711 logical ijnonbond 712 713 rmsd = sqrt(sum((xyz-nlist%hbrefgeo)**2))/dble(n) 714 715 if(rmsd.lt.1.d-6 .or. rmsd.gt. 0.3d0) then ! update list if first call or substantial move occured 716 717 nlist%nhb1=0 718 nlist%nhb2=0 719 do ix=1,topo%nathbAB 720 i=topo%hbatABl(1,ix) 721 j=topo%hbatABl(2,ix) 722 ij=j+i*(i-1)/2 723 rab=sqrab(ij) 724 if(rab.gt.hbthr1)cycle 725 ijnonbond=topo%bpair(ij).ne.1 726 do k=1,topo%nathbH 727 nh=topo%hbatHl(k) 728 inh=lin(i,nh) 729 jnh=lin(j,nh) 730 if(topo%bpair(inh).eq.1.and.ijnonbond)then ! exclude cases where A and B are bonded 731 nlist%nhb2=nlist%nhb2+1 732 nlist%hblist2(1,nlist%nhb2)=i 733 nlist%hblist2(2,nlist%nhb2)=j 734 nlist%hblist2(3,nlist%nhb2)=nh 735 elseif(topo%bpair(jnh).eq.1.and.ijnonbond)then ! exclude cases where A and B are bonded 736 nlist%nhb2=nlist%nhb2+1 737 nlist%hblist2(1,nlist%nhb2)=j 738 nlist%hblist2(2,nlist%nhb2)=i 739 nlist%hblist2(3,nlist%nhb2)=nh 740 elseif(rab+sqrab(inh)+sqrab(jnh).lt.hbthr2) then 741 nlist%nhb1=nlist%nhb1+1 742 nlist%hblist1(1,nlist%nhb1)=i 743 nlist%hblist1(2,nlist%nhb1)=j 744 nlist%hblist1(3,nlist%nhb1)=nh 745 endif 746 enddo 747 enddo 748 749 nlist%nxb =0 750 do ix=1,topo%natxbAB 751 i =topo%xbatABl(1,ix) 752 j =topo%xbatABl(2,ix) 753 ij=j+i*(i-1)/2 754 rab=sqrab(ij) 755 if(rab.gt.hbthr2)cycle 756 nlist%nxb=nlist%nxb+1 757 nlist%hblist3(1,nlist%nxb)=i 758 nlist%hblist3(2,nlist%nxb)=j 759 nlist%hblist3(3,nlist%nxb)=topo%xbatABl(3,ix) 760 enddo 761 762 nlist%hbrefgeo = xyz 763 764 endif ! else do nothing 765 766 end subroutine gfnff_hbset 767 768!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 769 770subroutine bond_hbset(n,at,xyz,sqrab,bond_hbn,bond_hbl,topo,hbthr1,hbthr2) 771use xtb_mctc_accuracy, only : wp 772 use xtb_gfnff_param 773 implicit none 774 type(TGFFTopology), intent(in) :: topo 775 integer,intent(in) :: n 776 integer,intent(in) :: at(n) 777 integer,intent(in) :: bond_hbn 778 integer,intent(out) :: bond_hbl(3,bond_hbn) 779 real(wp),intent(in) :: sqrab(n*(n+1)/2) 780 real(wp),intent(in) :: xyz(3,n) 781 real(wp),intent(in) :: hbthr1, hbthr2 782 783 integer i,j,k,nh,ia,ix,lin,ij,inh,jnh 784 integer bond_nr 785 real(wp) rab 786 logical ijnonbond 787 788 789 bond_nr=0 790 bond_hbl=0 791 do ix=1,topo%nathbAB 792 i=topo%hbatABl(1,ix) 793 j=topo%hbatABl(2,ix) 794 ij=j+i*(i-1)/2 795 rab=sqrab(ij) 796 if(rab.gt.hbthr1)cycle 797 ijnonbond=topo%bpair(ij).ne.1 798 do k=1,topo%nathbH 799 nh=topo%hbatHl(k) 800 inh=lin(i,nh) 801 jnh=lin(j,nh) 802 if(topo%bpair(inh).eq.1.and.ijnonbond)then ! exclude cases where A and B are bonded 803 bond_nr=bond_nr+1 804 bond_hbl(1,bond_nr)=i 805 bond_hbl(2,bond_nr)=j 806 bond_hbl(3,bond_nr)=nh 807 elseif(topo%bpair(jnh).eq.1.and.ijnonbond)then ! exclude cases where A and B are bonded 808 bond_nr=bond_nr+1 809 bond_hbl(1,bond_nr)=j 810 bond_hbl(2,bond_nr)=i 811 bond_hbl(3,bond_nr)=nh 812 endif 813 enddo 814 enddo 815 816end subroutine bond_hbset 817 818subroutine bond_hbset0(n,at,xyz,sqrab,bond_hbn,topo,hbthr1,hbthr2) 819use xtb_mctc_accuracy, only : wp 820 use xtb_gfnff_param 821 implicit none 822 type(TGFFTopology), intent(in) :: topo 823 integer,intent(in) :: n 824 integer,intent(in) :: at(n) 825 integer,intent(out) :: bond_hbn 826 real(wp),intent(in) :: sqrab(n*(n+1)/2) 827 real(wp),intent(in) :: xyz(3,n) 828 real(wp),intent(in) :: hbthr1, hbthr2 829 830 integer i,j,k,nh,ia,ix,lin,ij,inh,jnh 831 real(wp) rab 832 logical ijnonbond 833 834 835 bond_hbn=0 836 do ix=1,topo%nathbAB 837 i=topo%hbatABl(1,ix) 838 j=topo%hbatABl(2,ix) 839 ij=j+i*(i-1)/2 840 rab=sqrab(ij) 841 if(rab.gt.hbthr1)cycle 842 ijnonbond=topo%bpair(ij).ne.1 843 do k=1,topo%nathbH 844 nh=topo%hbatHl(k) 845 inh=lin(i,nh) 846 jnh=lin(j,nh) 847 if(topo%bpair(inh).eq.1.and.ijnonbond)then ! exclude cases where A and B are bonded 848 bond_hbn=bond_hbn+1 849 elseif(topo%bpair(jnh).eq.1.and.ijnonbond)then ! exclude cases where A and B are bonded 850 bond_hbn=bond_hbn+1 851 endif 852 enddo 853 enddo 854 855end subroutine bond_hbset0 856 857subroutine bond_hb_AHB_set(n,at,numbond,bond_hbn,bond_hbl,tot_AHB_nr,lin_AHB,topo) 858 use xtb_gfnff_param 859 implicit none 860 !Dummy 861 type(TGFFTopology), intent(inout) :: topo 862 integer,intent(in) :: n 863 integer,intent(in) :: numbond 864 integer,intent(in) :: at(n) 865 integer,intent(in) :: bond_hbn 866 integer,intent(in) :: bond_hbl(3,bond_hbn) 867 integer,intent(in) :: tot_AHB_nr 868 integer,intent(inout) :: lin_AHB(0:tot_AHB_nr) 869 !Stack 870 integer :: i,j 871 integer :: ii,jj 872 integer :: ia,ja 873 integer :: lin 874 integer :: hbH,hbA 875 integer :: Hat,Aat 876 integer :: Bat,atB 877 integer :: tot_count 878 integer :: AH_count 879 integer :: B_count 880 integer :: lin_diff 881 882 tot_count=0 883 AH_count=0 884 B_count=0 885 lin_diff=0 886 887 do i=1,numbond 888 ii=topo%blist(1,i) 889 jj=topo%blist(2,i) 890 ia=at(ii) 891 ja=at(jj) 892 if (ia.eq.1) then 893 hbH=ii 894 hbA=jj 895 else if (ja.eq.1) then 896 hbH=jj 897 hbA=ii 898 else 899 cycle 900 end if 901 if (at(hbA).eq.7.or.at(hbA).eq.8) then 902 do j = 1, bond_hbn 903 Bat = bond_hbl(2,j) 904 atB = at(Bat) 905 Aat = bond_hbl(1,j) 906 Hat = bond_hbl(3,j) 907 if (hbA.eq.Aat.and.hbH.eq.Hat) then 908 if (atB.eq.7.or.atB.eq.8) then 909 tot_count = tot_count + 1 910 lin_AHB(tot_count) = lin(hbA,hbH) 911 lin_diff=lin_AHB(tot_count)-lin_AHB(tot_count-1) 912 if (lin_diff.eq.0) then 913 B_count = B_count + 1 914 end if 915 !Next AH pair 916 if (lin_diff.ne.0) then 917 AH_count = AH_count + 1 918 topo%bond_hb_AH(1,AH_count) = hbA 919 topo%bond_hb_AH(2,AH_count) = hbH 920 !Reset B count 921 B_count = 1 922 end if 923 topo%bond_hb_Bn(AH_count) = B_count 924 topo%bond_hb_B(B_count,AH_count) = Bat 925 end if 926 else 927 cycle 928 end if 929 end do 930 end if 931 end do 932 933end subroutine bond_hb_AHB_set 934 935subroutine bond_hb_AHB_set1(n,at,numbond,bond_hbn,bond_hbl,tot_AHB_nr,lin_AHB,AH_count,bmax,topo) 936 use xtb_gfnff_param 937 implicit none 938 !Dummy 939 type(TGFFTopology), intent(inout) :: topo 940 integer,intent(in) :: n 941 integer,intent(in) :: numbond 942 integer,intent(in) :: at(n) 943 integer,intent(in) :: bond_hbn 944 integer,intent(in) :: bond_hbl(3,bond_hbn) 945 integer,intent(in) :: tot_AHB_nr 946 integer,intent(inout) :: lin_AHB(0:tot_AHB_nr) 947 integer,intent(out) :: AH_count 948 integer,intent(out) :: bmax 949 !Stack 950 integer :: i,j 951 integer :: ii,jj 952 integer :: ia,ja 953 integer :: lin 954 integer :: hbH,hbA 955 integer :: Hat,Aat 956 integer :: Bat,atB 957 integer :: tot_count 958 integer :: B_count 959 integer :: lin_diff 960 961 tot_count=0 962 AH_count=0 963 B_count=1 964 bmax=1 965 lin_diff=0 966 967 do i=1,numbond 968 ii=topo%blist(1,i) 969 jj=topo%blist(2,i) 970 ia=at(ii) 971 ja=at(jj) 972 if (ia.eq.1) then 973 hbH=ii 974 hbA=jj 975 else if (ja.eq.1) then 976 hbH=jj 977 hbA=ii 978 else 979 cycle 980 end if 981 if (at(hbA).eq.7.or.at(hbA).eq.8) then 982 do j = 1, bond_hbn 983 Bat = bond_hbl(2,j) 984 atB = at(Bat) 985 Aat = bond_hbl(1,j) 986 Hat = bond_hbl(3,j) 987 if (hbA.eq.Aat.and.hbH.eq.Hat) then 988 if (atB.eq.7.or.atB.eq.8) then 989 tot_count = tot_count + 1 990 lin_AHB(tot_count) = lin(hbA,hbH) 991 lin_diff=lin_AHB(tot_count)-lin_AHB(tot_count-1) 992 if (lin_diff.eq.0) B_count = B_count + 1 993 if (lin_diff.ne.0) then 994 AH_count = AH_count + 1 995 B_count = 1 996 end if 997 if (B_count.gt.bmax) bmax = B_count 998 end if 999 else 1000 cycle 1001 end if 1002 end do 1003 topo%nr_hb(i) = B_count 1004 end if 1005 end do 1006 1007end subroutine bond_hb_AHB_set1 1008 1009subroutine bond_hb_AHB_set0(n,at,numbond,bond_hbn,bond_hbl,tot_AHB_nr,topo) 1010 use xtb_gfnff_param 1011 implicit none 1012 !Dummy 1013 type(TGFFTopology), intent(in) :: topo 1014 integer,intent(in) :: n 1015 integer,intent(in) :: numbond 1016 integer,intent(in) :: at(n) 1017 integer,intent(in) :: bond_hbn 1018 integer,intent(in) :: bond_hbl(3,bond_hbn) 1019 integer,intent(out) :: tot_AHB_nr 1020 !Stack 1021 integer :: i,j 1022 integer :: ii,jj 1023 integer :: ia,ja 1024 integer :: hbH,hbA 1025 integer :: Hat,Aat 1026 integer :: Bat,atB 1027 1028 tot_AHB_nr=0 1029 1030 do i=1,numbond 1031 ii=topo%blist(1,i) 1032 jj=topo%blist(2,i) 1033 ia=at(ii) 1034 ja=at(jj) 1035 if (ia.eq.1) then 1036 hbH=ii 1037 hbA=jj 1038 else if (ja.eq.1) then 1039 hbH=jj 1040 hbA=ii 1041 else 1042 cycle 1043 end if 1044 if (at(hbA).eq.7.or.at(hbA).eq.8) then 1045 do j = 1, bond_hbn 1046 Bat = bond_hbl(2,j) 1047 atB = at(Bat) 1048 Aat = bond_hbl(1,j) 1049 Hat = bond_hbl(3,j) 1050 if (hbA.eq.Aat.and.hbH.eq.Hat) then 1051 if (atB.eq.7.or.atB.eq.8) then 1052 tot_AHB_nr = tot_AHB_nr + 1 1053 end if 1054 else 1055 cycle 1056 end if 1057 end do 1058 end if 1059 end do 1060 1061end subroutine bond_hb_AHB_set0 1062 1063!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1064 1065subroutine gfnff_hbset0(n,at,xyz,sqrab,topo,nhb1,nhb2,nxb,hbthr1,hbthr2) 1066use xtb_mctc_accuracy, only : wp 1067 use xtb_gfnff_param 1068 implicit none 1069 type(TGFFTopology), intent(in) :: topo 1070 integer, intent(out) :: nhb1 1071 integer, intent(out) :: nhb2 1072 integer, intent(out) :: nxb 1073 integer n 1074 integer at(n) 1075 real(wp) sqrab(n*(n+1)/2) 1076 real(wp) xyz(3,n) 1077 real(wp),intent(in) :: hbthr1, hbthr2 1078 1079 integer i,j,k,nh,ia,ix,lin,ij,inh,jnh 1080 logical ijnonbond 1081 real(wp) rab 1082 1083 nhb1=0 1084 nhb2=0 1085 do ix=1,topo%nathbAB 1086 i=topo%hbatABl(1,ix) 1087 j=topo%hbatABl(2,ix) 1088 ij=j+i*(i-1)/2 1089 rab=sqrab(ij) 1090 if(rab.gt.hbthr1)cycle 1091 ijnonbond=topo%bpair(ij).ne.1 1092 do k=1,topo%nathbH 1093 nh=topo%hbatHl(k) 1094 inh=lin(i,nh) 1095 jnh=lin(j,nh) 1096 if(topo%bpair(inh).eq.1.and.ijnonbond)then 1097 nhb2=nhb2+1 1098 elseif(topo%bpair(jnh).eq.1.and.ijnonbond)then 1099 nhb2=nhb2+1 1100 elseif(rab+sqrab(inh)+sqrab(jnh).lt.hbthr2) then 1101 nhb1=nhb1+1 1102 endif 1103 enddo 1104 enddo 1105 1106 nxb =0 1107 do ix=1,topo%natxbAB 1108 i =topo%xbatABl(1,ix) 1109 j =topo%xbatABl(2,ix) 1110 ij=j+i*(i-1)/2 1111 rab=sqrab(ij) 1112 if(rab.gt.hbthr2)cycle 1113 nxb=nxb+1 1114 enddo 1115 1116 end subroutine gfnff_hbset0 1117 1118!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1119! HB strength 1120!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1121 1122 subroutine hbonds(i,j,ci,cj,param,topo) 1123 use xtb_gfnff_param 1124 implicit none 1125 type(TGFFTopology), intent(in) :: topo 1126 type(TGFFData), intent(in) :: param 1127 integer i,j 1128 integer ati,atj 1129 real*8 ci(2),cj(2) 1130 ci(1)=topo%hbbas(i) 1131 cj(1)=topo%hbbas(j) 1132 ci(2)=topo%hbaci(i) 1133 cj(2)=topo%hbaci(j) 1134 end subroutine hbonds 1135 1136!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1137! ring analysis routine, don't touch ;-) 1138!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1139 1140subroutine getring36(n,at,nbin,a0_in,cout,irout) 1141 implicit none 1142 integer cout(10,20),irout(20) ! output: atomlist, ringsize, # of rings in irout(20) 1143 integer at(n) 1144 integer n,nbin(20,n),a0,i,nb(20,n),a0_in 1145 integer i1,i2,i3,i4,i5,i6,i7,i8,i9,i10 1146 integer n0,n1,n2,n3,n4,n5,n6,n7,n8,n9,n10 1147 integer a1,a2,a3,a4,a5,a6,a7,a8,a9,a10 1148 integer maxr 1149 parameter (maxr=500) 1150 integer list(n),m,mm,nn,c(10),cdum(10,maxr),iring 1151 integer adum1(0:n),adum2(0:n),kk,j,idum(maxr),same(maxr),k 1152 real*8 w(n),av,sd 1153 1154 cout = 0 1155 irout= 0 1156! if(n.le.2.or.at(a0_in).gt.10.or.nbin(19,a0_in).eq.1) return 1157 if(n.le.2. .or.nbin(19,a0_in).eq.1) return 1158 1159 nn=nbin(20,a0_in) 1160 1161 cdum=0 1162 kk=0 1163 do m=1,nn 1164! if(nb(m,a0_in).eq.1) cycle ! check (comment out) 1165 nb=nbin 1166 if(nb(m,a0_in).eq.1) cycle ! check (comment out) 1167 do i=1,n 1168 if(nb(20,i).eq.1)nb(20,i)=0 1169 enddo 1170 1171 do mm=1,nn 1172 w(mm)=dble(mm) 1173 list(mm)=mm 1174 enddo 1175 w(m) =0.0d0 1176 call ssort(nn,w,list) 1177 do mm=1,nn 1178 nb(mm,a0_in)=nbin(list(mm),a0_in) 1179 enddo 1180 1181 iring =0 1182 c =0 1183 1184 a0=a0_in 1185 n0=nb(20,a0) 1186 1187 do i1=1,n0 1188 a1=nb(i1,a0) 1189 if(a1.eq.a0) cycle 1190 n1=nb(20,a1) 1191 do i2=1,n1 1192 a2=nb(i2,a1) 1193 if(a2.eq.a1) cycle 1194 n2=nb(20,a2) 1195 do i3=1,n2 1196 a3=nb(i3,a2) 1197 n3=nb(20,a3) 1198 if(a3.eq.a2) cycle 1199 c(1)=a1 1200 c(2)=a2 1201 c(3)=a3 1202 if(a3.eq.a0.and.chkrng(n,3,c))then 1203 iring=3 1204 if(kk.eq.maxr) goto 99 1205 kk=kk+1 1206 cdum(1:iring,kk)=c(1:iring) 1207 idum(kk)=iring 1208 endif 1209 do i4=1,n3 1210 a4=nb(i4,a3) 1211 n4=nb(20,a4) 1212 if(a4.eq.a3) cycle 1213 c(4)=a4 1214 if(a4.eq.a0.and.chkrng(n,4,c))then 1215 iring=4 1216 if(kk.eq.maxr) goto 99 1217 kk=kk+1 1218 cdum(1:iring,kk)=c(1:iring) 1219 idum(kk)=iring 1220 endif 1221 do i5=1,n4 1222 a5=nb(i5,a4) 1223 n5=nb(20,a5) 1224 if(a5.eq.a4) cycle 1225 c(5)=a5 1226 if(a5.eq.a0.and.chkrng(n,5,c))then 1227 iring=5 1228 if(kk.eq.maxr) goto 99 1229 kk=kk+1 1230 cdum(1:iring,kk)=c(1:iring) 1231 idum(kk)=iring 1232 endif 1233 do i6=1,n5 1234 a6=nb(i6,a5) 1235 n6=nb(20,a6) 1236 if(a6.eq.a5) cycle 1237 c(6)=a6 1238 if(a6.eq.a0.and.chkrng(n,6,c))then 1239 iring=6 1240 if(kk.eq.maxr) goto 99 1241 kk=kk+1 1242 cdum(1:iring,kk)=c(1:iring) 1243 idum(kk)=iring 1244 endif 1245 enddo 1246 enddo 1247 enddo 1248 enddo 1249 enddo 1250 enddo 1251 1252 99 continue 1253 1254 enddo 1255 1256! compare 1257 same=0 1258 do i=1,kk 1259 do j=i+1,kk 1260 if(idum(i).ne.idum(j)) cycle ! different ring size 1261 if(same(j).eq.1 ) cycle ! already double 1262 adum1=0 1263 adum2=0 1264 do m=1,10 1265 i1=cdum(m,i) 1266 i2=cdum(m,j) 1267 adum1(i1)=1 1268 adum2(i2)=1 1269 enddo 1270 if(sum(abs(adum1-adum2)).ne.0) then 1271 same(j)=0 1272 else 1273 same(j)=1 1274 endif 1275 enddo 1276 enddo 1277 1278 m=0 1279 do i=1,kk 1280 if(same(i).eq.0) then 1281 m=m+1 1282 irout(m)=idum(i) ! number of atoms in ring m 1283 nn=idum(i) 1284 cout(1:nn,m)=cdum(1:nn,i) 1285 if(m.gt.19) then 1286 m=19 1287 goto 999 1288 endif 1289 endif 1290 enddo 1291999 irout(20)=m ! number of rings for this atom 1292 1293 return 1294 end subroutine getring36 1295 1296 subroutine ssort(n,edum,ind) 1297 implicit none 1298 integer n,ii,k,j,m,i,sc1 1299 real*8 edum(n),pp 1300 integer ind(n) 1301 1302 do 140 ii = 2, n 1303 i = ii - 1 1304 k = i 1305 pp= edum(i) 1306 do 120 j = ii, n 1307 if (edum(j) .gt. pp) go to 120 1308 k = j 1309 pp= edum(j) 1310 120 continue 1311 if (k .eq. i) go to 140 1312 edum(k) = edum(i) 1313 edum(i) = pp 1314 sc1=ind(i) 1315 ind(i)=ind(k) 1316 ind(k)=sc1 1317 140 continue 1318 1319 end subroutine ssort 1320 1321!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1322! neighbor only version of EEQ model 1323! included up to 1,4 interactions 1324!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1325 1326subroutine goedeckera(env,n,at,nb,pair,q,es,topo) 1327 use xtb_mctc_accuracy, only : wp 1328 use xtb_mctc_lapack, only : mctc_sytrf, mctc_sytrs 1329 implicit none 1330 character(len=*), parameter :: source = 'gfnff_ini2_goedeckera' 1331 type(TEnvironment), intent(inout) :: env 1332 type(TGFFTopology), intent(in) :: topo 1333 integer, intent(in) :: n ! number of atoms 1334 integer, intent(in) :: at(n) ! ordinal numbers 1335 integer, intent(in) :: nb(20,n) ! neighbors 1336 real(wp),intent(in) :: pair(n*(n+1)/2) 1337 real(wp),intent(out) :: q(n) ! output charges 1338 real(wp),intent(out) :: es ! ES energy 1339 1340! local variables 1341 logical :: exitRun 1342 integer :: m,i,j,k,l,ii,jj,kk 1343 integer :: ij,lj 1344 integer,allocatable :: ipiv(:) 1345 1346 real(wp) :: gammij,sief1,sief2 1347 real(wp) :: r2,r0 1348 real(wp) :: rij 1349 real(wp) :: tsqrt2pi,bohr 1350 real(wp) :: tmp 1351 real(wp),allocatable :: A (:,:) 1352 real(wp),allocatable :: x(:) 1353 1354! parameter 1355 parameter (tsqrt2pi = 0.797884560802866_wp) 1356 1357 m=n+topo%nfrag ! # atoms frag constrain 1358 allocate(A(m,m),x(m),ipiv(m)) 1359 1360! call prmati(6,pair,n,0,'pair') 1361 1362 A = 0 1363 1364! setup RHS 1365 do i=1,n 1366 x(i) =topo%chieeq(i) ! EN of atom 1367 A(i,i) =topo%gameeq(i)+tsqrt2pi/sqrt(topo%alpeeq(i)) 1368 enddo 1369 1370! setup A matrix 1371 do i=1,n 1372 do j=1,i-1 1373 ij = i*(i-1)/2+j 1374 rij=pair(ij) 1375 r2 = rij*rij 1376 gammij=1.d0/sqrt(topo%alpeeq(i)+topo%alpeeq(j)) ! squared above 1377 tmp = erf(gammij*rij)/rij 1378 A(j,i) = tmp 1379 A(i,j) = tmp 1380 enddo 1381 enddo 1382 1383! fragment charge constrain 1384 do i=1,topo%nfrag 1385 x(n+i)=topo%qfrag(i) 1386 do j=1,n 1387 if(topo%fraglist(j).eq.i) then 1388 A(n+i,j)=1 1389 A(j,n+i)=1 1390 endif 1391 enddo 1392 enddo 1393! call prmat(6,A,m,m,'A ini') 1394 1395 call mctc_sytrf(env, a, ipiv) 1396 call mctc_sytrs(env, a, x, ipiv) 1397 1398 call env%check(exitRun) 1399 if (exitRun) then 1400 call env%error('Solving linear equations failed', source) 1401 return 1402 end if 1403 1404 q(1:n) = x(1:n) 1405 1406 if(n.eq.1) q(1)=topo%qfrag(1) 1407 1408! energy 1409 es = 0.0_wp 1410 do i=1,n 1411 ii = i*(i-1)/2 1412 do j=1,i-1 1413 ij = ii+j 1414 rij=pair(ij) 1415 gammij=1.d0/sqrt(topo%alpeeq(i)+topo%alpeeq(j)) ! squared above 1416 tmp = erf(gammij*rij)/rij 1417 es = es + q(i)*q(j)*tmp/rij 1418 enddo 1419 es = es - q(i)* topo%chieeq(i) & 1420 & + q(i)*q(i)*0.5d0*(topo%gameeq(i)+tsqrt2pi/sqrt(topo%alpeeq(i))) 1421 enddo 1422 1423end subroutine goedeckera 1424 1425!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1426! condense charges to heavy atoms based on topology 1427!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1428 1429 subroutine qheavy(n,at,nb,q) 1430 implicit none 1431 integer,intent(in) :: n,nb(20,n),at(n) 1432 real*8,intent(inout) :: q(n) 1433 1434 integer i,j,k 1435 real*8 qtmp(n) 1436 1437 qtmp = q 1438 do i=1,n 1439 if(at(i).ne.1) cycle 1440 qtmp(i)=0 1441 do j=1,nb(20,i) 1442 k=nb(j,i) 1443 qtmp(k)=qtmp(k)+q(i)/dble(nb(20,i)) ! could be a bridging H 1444 enddo 1445 enddo 1446 1447 q = qtmp 1448 1449 end subroutine qheavy 1450 1451!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1452! determine number of cov. bonds between atoms 1453!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1454 1455 subroutine nbondmat(n,nb,pair) 1456 implicit none 1457! Dummy 1458 integer,intent(in) :: n 1459 integer,intent(in) :: nb(20,n) 1460 integer,intent(out) :: pair(n*(n+1)/2) 1461! Stack 1462 integer i,ni,newi,j,newatom,tag,d,i1,ni1,iii,ii,jj,k,lin 1463 integer,allocatable :: list(:,:),nlist(:,:),nnn(:),nn(:) 1464 logical da 1465 1466 allocate(nnn(n),nn(n),list(5*n,n),nlist(5*n,n)) 1467 1468 nn(1:n)=nb(20,1:n) 1469 1470 pair=0 1471 list=0 1472 do i=1,n 1473 ni=nn(i) 1474 list(1:ni,i)=nb(1:ni,i) 1475 enddo 1476 1477 nlist=list 1478 1479 pair=0 1480 do i=1,n 1481 do j=1,nb(20,i) 1482 k=nb(j,i) 1483 pair(lin(k,i))=1 1484 enddo 1485 enddo 1486 1487! one bond, tag=1 1488 tag=1 1489! call pairsbond(n,nn,list,pair,tag) 1490 1491! determine up to 3 bonds in between 1492 do d=1,2 1493 1494! loop over atoms 1495 do i=1,n 1496 ni=nn(i) 1497 newi=ni 1498! all neighbors of i 1499 do ii=1,ni 1500 i1=list(ii,i) 1501 ni1=nb(20,i1) 1502! all neighbors of neighbors of i 1503 do iii=1,ni1 1504 newatom=nb(iii,i1) 1505 da=.false. 1506 do j=1,newi 1507 if(newatom.eq.list(j,i))da=.true. 1508 enddo 1509 if(.not.da)then 1510 newi=newi+1 1511 nlist(newi,i)=newatom 1512 endif 1513 enddo 1514 enddo 1515 nnn(i)=newi 1516 enddo 1517 1518 list=nlist 1519 nn =nnn 1520 1521! one bond more 1522 tag=tag+1 1523 call pairsbond(n,nn,list,pair,tag) 1524 1525 enddo 1526 do i=1,n 1527 do j=1,i 1528 if(i.ne.j.and.pair(lin(j,i)).eq.0) pair(lin(j,i))=5 1529 enddo 1530 enddo 1531 1532 end subroutine nbondmat 1533 1534 subroutine pairsbond(n,nn,list,pair,tag) 1535 implicit none 1536 integer n,nn(n),list(5*n,n),tag 1537 integer i,j,k,ni,nj,ii,jj,ij 1538 integer pair(n*(n+1)/2) 1539 logical dai,daj 1540 1541 do i=1,n 1542 ni=nn(i) 1543 ij= i*(i-1)/2 1544 do j=1,i-1 1545 k = ij+j 1546 nj=nn(j) 1547 dai=.false. 1548 daj=.false. 1549 do ii=1,ni 1550 if(list(ii,i).eq.j)daj=.true. 1551 enddo 1552 do jj=1,nj 1553 if(list(jj,j).eq.i)dai=.true. 1554 enddo 1555 if(dai.and.daj.and.pair(k).eq.0) then 1556 pair(k)=tag 1557 endif 1558 enddo 1559 enddo 1560 1561 end subroutine pairsbond 1562 1563!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1564 1565 logical function pilist(ati) 1566 integer ati 1567 pilist=.false. 1568! if(ati.eq.5.or.ati.eq.6.or.ati.eq.7.or.ati.eq.8.or.ati.eq.9.or.ati.eq.16) pilist=.true. 1569 if(ati.eq.5.or.ati.eq.6.or.ati.eq.7.or.ati.eq.8.or.ati.eq.9.or.ati.eq.16.or.ati.eq.17) pilist=.true. 1570 end function pilist 1571 1572 logical function nofs(ati) 1573 integer ati 1574 nofs=.false. 1575! if(ati.eq.7.or.ati.eq.8.or.ati.eq.9.or.ati.eq.16) nofs=.true. 1576 if(ati.eq.7.or.ati.eq.8.or.ati.eq.9.or.ati.eq.16.or.ati.eq.17) nofs=.true. 1577 end function nofs 1578 1579 logical function xatom(ati) 1580 integer ati 1581 xatom=.false. 1582 if(ati.eq.17.or.ati.eq.35.or.ati.eq.53.or.& ! X in A-X...B 1583 & ati.eq.16.or.ati.eq.34.or.ati.eq.52.or.& 1584 & ati.eq.15.or.ati.eq.33.or.ati.eq.51) xatom=.true. 1585 end function xatom 1586 1587!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1588 1589 integer function ctype(n,at,nb,pi,a) 1590 integer n,a,at(n),nb(20,n),pi(n) 1591 integer i,no,j 1592 1593 ctype = 0 ! don't know 1594 1595 no=0 1596 do i=1,nb(20,a) 1597 j=nb(i,a) 1598 if(at(j).eq.8.and.pi(j).ne.0) no = no + 1 1599 enddo 1600 1601 if(no.eq.1.and.pi(a).ne.0) ctype = 1 ! a C=O carbon 1602 1603 end function ctype 1604 1605!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1606 1607 logical function alphaCO(n,at,hyb,nb,pi,a,b) 1608 integer n,a,b,at(n),hyb(n),nb(20,n),pi(n) 1609 integer i,j,no,nc 1610 1611 alphaCO = .false. 1612 if(pi(a).ne.0.and.hyb(b).eq.3.and.at(a).eq.6.and.at(b).eq.6) then 1613 no = 0 1614 do i=1,nb(20,a) 1615 j=nb(i,a) 1616 if(at(j).eq. 8.and.pi(j).ne.0.and.nb(20,j).eq.1) no = no +1 ! a pi =O on the C? 1617 enddo 1618 if(no.eq.1) then 1619 alphaCO = .true. 1620 return 1621 endif 1622 endif 1623 if(pi(b).ne.0.and.hyb(a).eq.3.and.at(b).eq.6.and.at(a).eq.6) then 1624 no = 0 1625 do i=1,nb(20,b) 1626 j=nb(i,b) 1627 if(at(j).eq. 8.and.pi(j).ne.0.and.nb(20,j).eq.1) no = no +1 ! a pi =O on the C? 1628 enddo 1629 if(no.eq.1) then 1630 alphaCO = .true. 1631 return 1632 endif 1633 endif 1634 1635 end function alphaCO 1636 1637!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1638 1639 logical function amide(n,at,hyb,nb,pi,a) 1640 integer n,a,at(n),hyb(n),nb(20,n),pi(n) 1641 integer i,j,no,nc,ic 1642 1643 amide = .false. ! don't know 1644 if(pi(a).eq.0.or.hyb(a).ne.3.or.at(a).ne.7) return 1645 1646 nc=0 1647 no=0 1648 do i=1,nb(20,a) 1649 j=nb(i,a) 1650 if(at(j).eq.6.and.pi(j).ne.0) then ! a pi C on N? 1651 nc = nc + 1 1652 ic = j 1653 endif 1654 enddo 1655 1656 if(nc.eq.1)then 1657 do i=1,nb(20,ic) 1658 j=nb(i,ic) 1659 if(at(j).eq. 8.and.pi(j).ne.0.and.nb(20,j).eq.1) no = no +1 ! a pi =O on the C? 1660! if(at(j).eq.16.and.pi(j).ne.0.and.nb(20,j).eq.1) no = no +1 ! a pi =S on the C? 1661 enddo 1662 endif 1663 1664 if( no .eq. 1 ) amide = .true. 1665 1666 end function amide 1667 1668!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1669 1670 logical function amideH(n,at,hyb,nb,pi,a) 1671 integer n,a,at(n),hyb(n),nb(20,n),pi(n) 1672 integer i,j,nc,nn 1673 !logical amide 1674 1675 amideH = .false. ! don't know 1676 if(nb(20,a).ne.1 ) return 1677 nn=nb(1,a) ! the N 1678 if(.not.amide(n,at,hyb,nb,pi,nn)) return 1679 1680 nc=0 1681 do i=1,nb(20,nn) 1682 j=nb(i,nn) 1683 if(at(j).eq.6.and.hyb(j).eq.3) then ! a sp3 C on N? 1684 nc = nc + 1 1685 endif 1686 enddo 1687 1688 if( nc .eq. 1 ) amideH = .true. 1689 1690 end function amideH 1691 1692end module xtb_gfnff_ini2 1693