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 Licen 16! along with xtb. If not, see <https://www.gnu.org/licenses/>. 17 18 19! ..................................................................... 20! 21! Calculates the pointgroup, the nuclear exchange table and the 22! transfomation matrices 23! 24! ..................................................................... 25 26Subroutine symtrans(sflies,thrsym,xyz,natoms,ict,trans,ntrans) 27 use xtb_mctc_accuracy, only : wp 28 29 ! ..................................................................... 30 ! 31 ! Input 32 ! ----- 33 ! 34 ! xyz : cartesian coordinates of the atoms 35 ! 36 ! natoms : number of the atoms 37 ! 38 ! thrsym : threshhold 39 ! 40 ! 41 ! Output 42 ! ------ 43 ! 44 ! sflies : Schoenflies symbol 45 ! 46 ! ict : nuclear exchange table 47 ! 48 ! ntrans : number of transformations 49 ! 50 ! trans : transformation matrices 51 ! 52 ! ..................................................................... 53 54 implicit none 55 56 integer ndi14 57 parameter (ndi14=120) 58 59 logical show 60 61 character sflies*(*),csf*1,cla*1 62 63 integer natoms,nn,nprisy,iback,ngen,ntrans,mi(ndi14),& 64 & mj(ndi14),natomsred,nhessred,& 65 & invt(ndi14),ict(natoms,ndi14),ierr,i,j,k,l,& 66 & jkseq(2,natoms),ictnr(natoms,ndi14),ictt,& 67 & lieff(natoms) 68 69 real(wp) xyz(3,natoms),gen(9,3),trans(9,ndi14),thrsym 70 71 72 nprisy=-1 73 show=.false. 74 !sg 75 ict = 0 76 77 ! ... Decompose input Schoenflies symbol 78 call grpsmb(sflies,csf,cla,nn,nprisy) 79 80 ! ... Check validity of Schoenflies symbol 81 if (.not.(csf.eq.'s'.and.cla.eq.' '& 82 & .or. csf.eq.'c'.and.(cla.eq.'h'.or.cla.eq.'v'.or.cla.eq.' ')& 83 & .or. csf.eq.'d'.and.(cla.eq.'h'.or.cla.eq.'d'.or.cla.eq.' ')& 84 & .or. csf.eq.'t'.and.(cla.eq.'h'.or.cla.eq.'d'.or.cla.eq.' ')& 85 & .or. csf.eq.'o'.and.(cla.eq.'h'.or.cla.eq.' ')& 86 & .or. csf.eq.'i'.and.(cla.eq.'h'.or.cla.eq.' '))) then 87 ! write (*,'(a)') '<symmetry> INVALID SCHOENFLIES SYMBOL ! ' 88 ! write (*,*) '*',sflies,'*',csf,cla,nn 89 call flush(6) 90 ! error STOP 91 end if 92 93 ! ... Form generators of group corresponding to Schoenflies symbol 94 call getgen(csf,nn,cla,gen,ngen,iback,nprisy) 95 if(iback.ne.0) then 96 write (*,'(a)') '<symmetry> : Error in getgen' 97 call flush(6) 98 error STOP 99 end if 100 101 ! ... Generate transformation matrices 'trans' of the symmetry group 102 ! ... Find inverse group operations 'invt' and memorize group 103 ! ... generation path 'mi','mj' for subsequent construction 104 ! ... of representations 105 call group(gen,ngen,trans,ntrans,ndi14,thrsym,invt,mi,mj,iback,& 106 & nprisy) 107 108 ! ... 109 ! ... determines the nuclear exchange group 110 ! ... ict(i,ig) is the nucleus obtained on nucleus i 111 ! ... by means of the ig's symmetry operation 112 ! ... symmetries specified as 3*3 matrices on array trans 113 ! ... ntrans = group order 114 ! ... natoms = number of atoms 115 ! ... xyz = nuclear coordinates 116 ! ... 117 call nucex(ict,natoms,trans,ntrans,thrsym,xyz,natoms,ierr,show) 118 119 ! ... delete symmetry redundant atoms in ict 120 ictnr=ict 121 do i=1,natoms 122 do k=2,ntrans 123 ictt=ict(i,k) 124 do l=1,(k-1) 125 if (ictt.eq.ict(i,l)) ictnr(i,k)=0 126 end do 127 end do 128 end do 129 130 ! ... get list of sym.equivalent atoms jkseq(.,i), i=1,natoms 131 do i=1,natoms 132 jkseq(1,i)=0 133 jkseq(2,i)=0 134 do j=1,(i-1) 135 do k=2,ntrans 136 if (i.eq.ict(j,k).and.jkseq(1,i).eq.0) then 137 jkseq(1,i)=j 138 jkseq(2,i)=k 139 exit 140 end if 141 end do 142 end do 143 end do 144 natomsred=0 145 do i=1,natoms 146 if (jkseq(1,i).eq.0) then 147 natomsred=natomsred+1 148 lieff(natomsred)=i 149 end if 150 end do 151 nhessred=3*natomsred 152 !c write (*,*) 'natomsred=',natomsred 153 !c write (*,*) 'lieff(.)=',(lieff(i),i=1,natomsred) 154 155 !c write(*,*) 'sflies',sflies,'/' 156 !c write(*,*) 'ict ntrans=',ntrans 157 !c do i=1,natoms 158 !c write(*,*)(ict(i,j),j=1,ntrans) 159 !c end do 160 !c write(*,*) 'trans' 161 !c do j=1,ntrans 162 !c write(*,'(9f7.4)')(trans(i,j),i=1,9) 163 !c end do 164 !c write(*,*) 'ictnr',ntrans 165 !c do i=1,natoms 166 !c write(*,*)(ictnr(i,j),j=1,ntrans) 167 !c end do 168 !c write(*,*) 'jkseq' 169 !c do i=1,natoms 170 !c write(*,*)(jkseq(j,i),j=1,2) 171 !c end do 172 173End subroutine 174 175 176! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 177! $Id: grpsmb.f,v 1.21998/11/0315:21:49 klaus Exp $ 178! $Log: grpsmb.f,v $ 179! Revision 1.2 1998/11/0315:21:49 klaus 180! change getcor into allocate (F90) 181! 182! Revision 1.1 1992/09/1013:41:20 tomjones 183! Initial revision 184! 185! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 186!-------------------------------------------------------------------- 187 subroutine grpsmb(sflies,csf,cla,nn,nprisy) 188 use xtb_mctc_accuracy, only : wp 189 implicit real(wp)(a-h,o-z) 190 character csf,cla,zn(0:9),sflies*4 191 data zn /'0','1','2','3','4','5','6','7','8','9'/ 192! 193 nn=0 194! nn is the symmetry number of the z-axis 195 csf=sflies(1:1) 196 do i=1,9 197 if (sflies(2:2).eq.zn(i)) nn=i 198 end do 199 if (nn.eq.0) then 200 cla=sflies(2:2) 201! if (nprisy.ge.-1) write (6,'(/, 202! 1 '' symmetry group of the molecule : '',2a1)') csf,cla 203 if (csf.eq.'c'.and.(cla.eq.'s'.or.(nn.eq.0.and.cla.eq.'h')))& 204 & then 205 nn=1 206 cla='h' 207 elseif (csf.eq.'c' .and. cla.eq.'i') then 208 nn=2 209 cla=' ' 210 csf='s' 211 endif 212 else 213 np=-1 214 do i=0,9 215 if (sflies(3:3).eq.zn(i)) np=i 216 end do 217 if (np.ge.0) then 218 nn=10*nn+np 219 cla=sflies(4:4) 220! if (nprisy.ge.-1) 221! 1 write (6,'(/,'' symmetry group of the molecule : '', 222! 2 a1,i2,a1)') csf,nn,cla 223 else 224 cla=sflies(3:3) 225! if (nprisy.ge.-1) 226! 1 write (6,'(/,'' symmetry group of the molecule : '', 227! 2 a1,i1,a1)') csf,nn,cla 228 endif 229 endif 230 return 231 end subroutine 232! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 233! $Id: group.f,v 1.31998/11/0315:21:36 klaus Exp $ 234! $Log: group.f,v $ 235! Revision 1.3 1998/11/0315:21:36 klaus 236! change getcor into allocate (F90) 237! 238! Revision 1.2 1997/10/0113:07:02 marcok 239! The hp780 ld complained about not initialized variables. Done now. 240! 241! Revision 1.1 1992/09/1013:41:12 tomjones 242! Initial revision 243! 244! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 245!-------------------------------------------------------------------- 246subroutine group(gen,ngen,u,nt,max14,thrsym,inv,ipath,jpath,ierr,& 247 & nprisy) 248 use xtb_mctc_accuracy, only : wp 249 implicit real(wp)(a-h,o-z) 250 dimension u(9,max14),gen(9,ngen),inv(max14),scr(9) 251 dimension ipath(max14),jpath(max14) 252 ! generates 3x3 matrices forming group 253 ! from ngen generators on gen 254 ! output: 255 ! u = 3x3 matrices 256 ! nt = number of symmetry operations 257 ! max14 = maximum order of group 258 ! inv containes index of invers operation 259 it=0 260 a0=0 261 a1=1 262 ierr=0 263 ! set neutral element 264 ! (every point group has the neutral element as first operation) 265 nt = 1 266 do i=1,9 267 u(i,nt) = a0 268 end do 269 u(1,nt) = a1 270 u(5,nt) = a1 271 u(9,nt) = a1 272 ipath(1)=0 273 jpath(1)=0 274 if (ngen.eq.0) goto 700 275 do 500 igen=1,ngen 276 ! loop over generators 277 ncoset = 1 278 nto = nt 279 do i=1,9 280 scr(i) = gen(i,igen) 281 end do 282 ! check if generator is already contained in (sub)group 283 ! of order nto : 284 do j=1,nt 285 if(dif(u(1,j),scr(1),9).lt.thrsym) go to 500 286 end do 287 ipath(nt+1)=-igen 288 jpath(nt+1)=0 289 ! for new group element form coset : 290 245 do 250 it=1,nto 291 nt = nt + 1 292 if (nt.gt.max14) goto 900 293 if (it.eq.1) goto 250 294 ipath(nt)=nto*ncoset+1 295 jpath(nt)=it 296 250 call mult3 (u(1,nto*ncoset+it),scr,u(1,it),3) 297 ncoset = ncoset + 1 298 ! check : do subgroup plus cosets form a group ? 299 do 330 it = 1,nt 300 do 320 jcoset = 1,ncoset-1 301 call mult3 (scr,u(1,it),u(1,1+nto*jcoset),3) 302 do kt = 1,nt 303 if (dif(u(1,kt),scr(1),9).lt.thrsym) goto 320 304 end do 305 ! new symmetry operation found 306 ipath(nt+1)=it 307 jpath(nt+1)=1+nto*jcoset 308 goto 245 309 320 continue 310 330 continue 311 ! cosets now form group - take next generator 312 500 continue 313 if (nprisy.ge.2) then 314 write(6,'(/,i5,a)') nt,' symmetry operations found :' 315 elseif (nprisy.ge.0) then 316 write(6,'(/,i5,a)') nt,' symmetry operations found' 317 endif 318 ! 319 ! find inverse operators 320 inv(1) = 1 321 if (nt.eq.1) goto 700 322 do 660 i=2,nt 323 inv(i) = 0 324 do 640 j=2,nt 325 call mult3 (scr,u(1,i),u(1,j),3) 326 s = dif (scr,u(1,1),9) 327 if (s.gt.thrsym) go to 640 328 inv(i) = j 329 go to 660 330 640 continue 331 go to 950 332 660 continue 333 700 continue 334 ! 335 if (nprisy.ge.2) then 336 write(6,750) 337 750 format(/,10x,'xx',6x,'yx',6x,'zx',6x,'xy',6x,'yy',6x,& 338 & 'zy',6x,'xz',6x,'yz',6x,'zz') 339 do 5 it=1,nt 340 5 write(6,'(i4,1x,9f8.4)') it,(u(i,it),i=1,9) 341 write(6,*) 342 write(6,'(/,'' group element generation record'',& 343 & '' (negative numbers correspond to generators) :'')') 344 write (6,'(6(2x,i3,''='',i3,''*'',i3))')& 345 & (it,ipath(it),jpath(it),it=1,nt) 346 end if 347 ! 348 return 349 ! 350 900 write(6,'(a,i4)') ' size of group larger than parameter ndi14 =',& 351 & max14 352 ierr=1 353 return 354 ! 355 950 write(6,*) ' inversion problem ' 356 ierr=1 357 return 358end subroutine 359! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 360! $Id: nucex.f,v 1.42000/10/2013:51:18 haettig Exp $ 361! $Log: nucex.f,v $ 362! Revision 1.4 2000/10/2013:51:18 haettig 363! replaced fortran STOPs by call to quit routine. 364! 365! Revision 1.3 1998/11/0315:22:26 klaus 366! change getcor into allocate (F90) 367! 368! Revision 1.2 1995/04/0517:14:37 holger 369! Source formatting. 370! 371! Revision 1.1 1992/09/10 13:42:02 tomjones 372! Initial revision 373! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 374subroutine nucex(ict,max10,trans,ntrans,thrsym,xyz,natoms,& 375 & ierr,show) 376 use xtb_mctc_accuracy, only : wp 377! 378! determines the nuclear exchange group 379! ict(i,ig) is the nucleus obtained on nucleus i 380! by means of the ig's symmetry operation 381! symmetries specified as 3*3 matrices on array trans 382! ntrans = group order 383! natoms = number of atoms 384! xyz = nuclear coordinates 385! 386 implicit real(wp)(a-h,o-z) 387 dimension ict(max10,ntrans),trans(9,ntrans),scr(3),xyz(3,natoms) 388 logical show 389! 390! set neutral operation 391! 392 ierr=0 393 imotz=ierr 394!sg ierr=0 395 do 120 ic=1,natoms 396 ict(ic,1) = ic 397 120 continue 398 if (ntrans.eq.1) goto 400 399 do 300 ic=1,natoms 400 do 250 it=2,ntrans 401 ict(ic,it) = 0 402 call mult3 (scr,trans(1,it),xyz(1,ic),1) 403 do 200 jc=1,natoms 404 if (dif(scr,xyz(1,jc),3).gt.thrsym) goto 200 405 ict(ic,it) = jc 406 goto 250 407 200 continue 408 goto 500 409 250 continue 410 300 continue 411 400 continue 412!sg 413! open(unit=66,file='trans',form='unformatted') 414! write(66)ntrans 415! do it=1,ntrans 416! write(66) (trans(k,it),k=1,9) 417! enddo 418! do iatom=1,natoms 419! write(66) (ict(iatom,it),it=1,ntrans) 420! enddo 421! close(66) 422!sg 423 if (show) then 424 write(6,'(/,10x,a,/)') 'nuclear exchange table :' 425 write(6,'(1x,a2,2x,24i3)') 'it',(it,it=1,ntrans) 426 write(6,*) 427 do 30 iatom=1,natoms 428 write(6,'(5x,24i3)') (ict(iatom,it),it=1,ntrans) 429 30 continue 430 endif 431 return 432 500 ierr=1 433 if (imotz.eq.-1) return 434 write(6,*) ' error in nuclear exchange group determination' 435 write(6,*) ' coordinates ' 436 do 600 ic=1,natoms 437 write(6,'(1x,i5,5x,3f14.8)') ic,(xyz(k,ic),k=1,3) 438 600 continue 439 write(6,*) ' transformation matrices ' 440 do 700 it=1,ntrans 441 write(6,'(1x,i5)') it 442 write(6,'(1x,9f8.4)') (trans(k,it),k=1,9) 443 700 continue 444 error STOP 'fatal error in nucex.' 445 end subroutine 446! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 447! $Id: getgen.f,v 1.11992/09/1013:40:37 tomjones Exp $ 448! $Log: getgen.f,v $ 449! Revision 1.1 1992/09/1013:40:37 tomjones 450! Initial revision 451! 452! RCS cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 453!-------------------------------------------------------------------- 454 subroutine getgen (csf,nn,cla,gen,ngen,ierr,nprisy) 455 use xtb_mctc_accuracy, only : wp 456 implicit real(wp)(a-h,o-y) 457 implicit character (z) 458 dimension gen(9,3),unit(9) 459 character csf,cla 460 data zc,zd,zs,zt,zo,zi,zv,zh/'c','d','s','t','o','i','v','h'/ 461 a0=0 462 a1=1 463 a2=2 464 a3=3 465 a4=4 466 a5=5 467 ierr=0 468 pi=a4*atan(a1) 469 twopi=pi+pi 470! generate unit matrix 471 do 10 i=2,8 472 10 unit(i) = a0 473 unit(1) = a1 474 unit(5) = a1 475 unit(9) = a1 476 ngen=0 477! write(6,*) ' this program sets generators for point groups' 478! write(6,*) ' specified by schoenflies symbols, e.g. d2d' 479! write(6,*) ' please use s1 for cs and s2 for ci' 480! write(6,*) ' main axis is always the z - axis' 481! write(6,*) ' secondary axis is always x, e.g. for d3' 482! write(6,*) ' sigma(v) plane is always xz, e.g. for c3v' 483! write(6,*) ' input schoenflies symbol: c2v or td etc' 484 if (nprisy.ge.0) write(6,20) 485 20 format(/,' the group has the following generators :') 486! note that in ih the mirror plane is perpendicular to the c2-axis 487! but not to the c5- or c3-axis 488 if (csf.eq.zi.and.cla.eq.zh) cla=zv 489! 490! set generators for input group 491! 492! cn(z) 493! 494 n=0 495 if (nn.gt.0.and.(csf.eq.zc.or.csf.eq.zd).and.cla.ne.zd) n=nn 496 if (csf.eq.zt.and.cla.ne.zd) n=2 497 if (csf.eq.zo) n=4 498 if(csf.eq.zi) n=5 499 if (n.le.0) go to 250 500 if (n.le.9 .and. nprisy.ge.0) write(6,30) n 501 if (n.gt.9 .and. nprisy.ge.0) write(6,31) n 502 30 format(3x,'c',i1,'(z)') 503 31 format(3x,'c',i2,'(z)') 504 ngen = ngen +1 505 if (ngen.gt.3) go to 1000 506 do 200 i=1,9 507 200 gen(i,ngen) = unit(i) 508 an = n 509 alpha = twopi/an 510 x = sin(alpha) 511 y = cos(alpha) 512 gen(1,ngen) = y 513 gen(5,ngen) = y 514 gen(4,ngen) = -x 515 gen(2,ngen) = x 516 250 continue 517! 518! sn(z) 519! 520 n=0 521 if (csf.eq.zs) n=nn 522 if (csf.eq.cla) n=2*nn 523 if (csf.eq.zt.and.cla.eq.zd) n=4 524 if (n.le.0) go to 400 525 if (n.le.9 .and. nprisy.ge.0) write(6,300) n 526 if (n.gt.9 .and. nprisy.ge.0) write(6,301) n 527 300 format(3x,'s',i1,'(z)') 528 301 format(3x,'s',i2,'(z)') 529 ngen = ngen +1 530 if (ngen.gt.3) go to 1000 531 do 350 i=1,9 532 350 gen(i,ngen) =-unit(i) 533 an = n 534 alpha = twopi/an 535 x = sin(alpha) 536 y = cos(alpha) 537 gen(1,ngen) = y 538 gen(5,ngen) = y 539 gen(4,ngen) = -x 540 gen(2,ngen) = x 541 400 continue 542! 543! c3(1,1,1) for t,td,th,o,oh 544! 545 if (csf.ne.zt.and.csf.ne.zo) go to 450 546 if (nprisy.ge.0) write(6,410) 547 410 format(3x,'c3(1,1,1)') 548 ngen = ngen +1 549 if (ngen.gt.3) go to 1000 550 do 420 i=1,9 551 420 gen(i,ngen) = a0 552 gen(2,ngen) = a1 553 gen(6,ngen) = a1 554 gen(7,ngen) = a1 555 450 continue 556! 557! c3 for i or ih 558! 559 if (csf.ne.zi) go to 500 560 if (nprisy.ge.0) write(6,480) 561 480 format(3x,'c3 for i or ih') 562 ngen = ngen + 1 563 if (ngen.gt.3) go to 1000 564 sqr3 = sqrt(a3) 565 cosd = sin(a2*pi/a5)/((a1-cos(a2*pi/a5))*sqr3) 566 sind = sqrt(a1-cosd**2) 567 gen(1,ngen) = a1 - a3*cosd**2/a2 568 gen(2,ngen) = sqr3*cosd/a2 569 gen(3,ngen) = a3*sind*cosd/a2 570 gen(4,ngen) = -gen(2,ngen) 571 gen(5,ngen) = -a1/a2 572 gen(6,ngen) = sqr3*sind/a2 573 gen(7,ngen) = gen(3,ngen) 574 gen(8,ngen) = -gen(6,ngen) 575 gen(9,ngen) = a1-a3*sind**2/a2 576 500 continue 577! 578! c2(x) 579! 580 if(csf.ne.zd) go to 600 581 if (nprisy.ge.0) write(6,520) 582 520 format(3x,'c2(x)') 583 i=1 584 ngen = ngen +1 585 if (ngen.gt.3) go to 1000 586 do 550 j=1,9 587 550 gen(j,ngen) = -unit(j) 588 gen(4*i-3,ngen) = a1 589 600 continue 590! 591! mirror plane sigma(xz) or sigma(xy) 592! 593 if (cla.ne.zv.and.cla.ne.zh) go to 750 594 if (cla.eq.zv) then 595 if (nprisy.ge.0) write(6,610) 596 610 format(3x,'mirror plane sigma(xz)') 597 i=2 598 else 599 if (nprisy.ge.0) write(6,620) 600 620 format(3x,'mirror plane sigma(xy)') 601 i=3 602 endif 603 ngen = ngen + 1 604 if (ngen.gt.3) go to 1000 605 do 700 j=1,9 606 700 gen(j,ngen) = unit(j) 607 gen(4*i-3,ngen) = -a1 608 750 return 609 1000 write(6,*) ' more than ',3,' generators' 610 ierr=1 611 end subroutine 612 613 614 subroutine symmetry(natoms,nat3,xyz0,h,eig,totsym,nvar) 615 use xtb_mctc_accuracy, only : wp 616 implicit none 617 integer natoms,nat3,totsym(nat3),nvar 618 real(wp) h(nat3,nat3),xyz0(3,natoms),eig(nat3) 619 real(wp) :: xyz(3,natoms) 620 621 real(wp),allocatable :: trans(:,:) 622 real(wp) thrsym 623 integer ntrans,i,j,nat,it,k,m 624 logical lok 625 integer :: ifile 626 627 call open_binary(ifile,'trans','r') 628 read(ifile)ntrans 629 allocate(trans(9,ntrans)) 630 do it=1,ntrans 631 read(ifile) (trans(i,it),i=1,9) 632 enddo 633 call close_file(ifile) 634 635 thrsym=1.d-6*natoms 636 totsym = 0 637 nvar = 0 638 do i=1,nat3 639 xyz = xyz0 640 m = 0 641 do j=1,natoms 642 do k=1,3 643 m = m + 1 644 xyz(k,j)=xyz(k,j)+h(m,i)*0.005 645 enddo 646 enddo 647 call checksym(thrsym,natoms,ntrans,trans,xyz,lok) 648 if(lok.and.abs(eig(i)).gt.1.d-8)then 649 totsym(i)=1 650 nvar = nvar + 1 651 endif 652 enddo 653 write(*,*)'Number of symmetry operations ',ntrans 654 write(*,*)'Nvar, symmetry restricted ',nvar 655 656 end subroutine 657 658!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 659 660 subroutine grdsym(grd,natoms) 661 use xtb_mctc_accuracy, only : wp 662 use xtb_symparam 663!,ntrans,ict,trans) 664 665 implicit real(wp) (a-h,o-z) 666 dimension grd(3,natoms) 667! dimension ict(natoms,ntrans),trans(9,ntrans) 668 data garnix/1.d-14/, zero/0.d0/, one/1.d0/ 669 real(wp) :: det(3,natoms),scr(3) 670 671! ----- symmetrize cartesian gradient/displ vector ----- 672! using transformation matrices and nuclear exchange table ict 673! ict = transformation table "atoms versus symmetry operations" 674 675 if (ntrans.gt.1) then 676 do 200 iat = 1,natoms 677 det(1,iat) = zero 678 det(2,iat) = zero 679 det(3,iat) = zero 680 200 continue 681 do 300 iat=1,natoms 682 do 400 itrans=1,ntrans 683 call mult3(scr,trans(1,itrans),grd(1,iat),1) 684 isymat=ict(iat,itrans) 685 det(1,isymat)=det(1,isymat) + scr(1) 686 det(2,isymat)=det(2,isymat) + scr(2) 687 det(3,isymat)=det(3,isymat) + scr(3) 688 400 continue 689 300 continue 690 ogroup=one/dble(ntrans) 691 do 500 iat=1,natoms 692 grd(1,iat)=det(1,iat)*ogroup 693 grd(2,iat)=det(2,iat)*ogroup 694 grd(3,iat)=det(3,iat)*ogroup 695 500 continue 696 end if 697 698! brush away numerical zeros 699 do 575 iat=1,natoms 700 do 575 j=1,3 701 if (dabs(grd(j,iat)).lt.garnix) grd(j,iat) = zero 702 575 continue 703 704 return 705 end subroutine 706 707!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 708 709 subroutine checksym(thrsym,natoms,ntrans,trans,xyz,same) 710 use xtb_mctc_accuracy, only : wp 711 implicit none 712 real(wp) xyz(3,natoms),trans(9,ntrans) 713 integer natoms,ntrans 714 integer ic,it,jc 715 logical same 716 real(wp) scr(3),thrsym 717 real(wp) dif 718 719 same=.true. 720 do 300 ic=1,natoms 721 do 250 it=2,ntrans 722 call mult3 (scr,trans(1,it),xyz(1,ic),1) 723 do 200 jc=1,natoms 724 if (dif(scr,xyz(1,jc),3).gt.thrsym) goto 200 725 goto 250 726 200 continue 727 same = .false. 728 return 729 250 continue 730 300 continue 731 732 end subroutine 733 734! --------------------------------------------------------------------- 735 736 function dif(a,b,n) 737 use xtb_mctc_accuracy, only : wp 738 implicit real(wp) (a-h,o-z) 739 dimension a(n),b(n) 740 s = 0.d0 741 do 10 k=1,n 742 10 s = s + (b(k)-a(k))*(b(k)-a(k)) 743 dif = dsqrt(s) 744 return 745 end function 746 747! --------------------------------------------------------------------- 748 749 subroutine mult3 (a,b,c,n) 750 use xtb_mctc_accuracy, only : wp 751 implicit real(wp) (a-h,o-z) 752! --------------------------------------------------------------------- 753! multiplication of two matrices b,c 754! a = b*c dim: a(3,n),b(3,3),c(3,n) 755! --------------------------------------------------------------------- 756 dimension a(3,n),b(3,3),c(3,n) 757 do 30 j=1,n 758 do 20 i=1,3 759 s = 0.d0 760 do 10 k=1,3 761 10 s = s + b(i,k)*c(k,j) 762 20 a(i,j) = s 763 30 continue 764 765 return 766 end subroutine 767 768