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