#ifndef MY_TINY #define MY_TINY 1E-5_q #endif !********************* Symmetry package ******************************** ! * ! Author: Juergen Furthmueller * ! Max-Planck-Institut fuer Metallforschung * ! Heissenbergstr. 1 * ! D 7000 Stuttgart 80 * ! * ! Tel.: Germany-711 6860 533 * ! * !*********************************************************************** ! * SUBROUTINE SETGRP(S,SYMOP,GTRANS,NROT,NROTK,TAU,LATTYP,NOP, & & TAUROT,NSP,NSX,NA,NAX,INDEX,WRK,PRCHAN) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine SETGRP is the primary user-interface of this symmetry * ! package. The function of SETGRP is to set up all possible * ! space group operators (INTEGER rotation matrices and nontrivial * ! translations given in crystal coordinates) of a lattice with * ! some arbitrary basis (atomic arrangement). This will be done by * ! first setting up the point group operators of the pure bravais * ! lattice without basis (empty supercell) and checking which of the * ! operations can reproduce the lattice with basis (to be done in * ! routine GETGRP - for further details see comments on GETGRP). * ! * ! * ! Input parameters: * ! ----------------- * ! * ! TAU(NAX,3,NSX) contains the crystal coordinates of all atomic * ! positions within the supercell. Warning: TAU will be * ! modified on output! TAU will contain the same set of * ! atomic position as on input, but the order of storage * ! of the positions will be changed in some way described * ! in routine LATORD of the lattice vector package and * ! periodic boundary conditions are imposed on the atomic * ! coordinates to shift all atoms into a box defined by * ! the corner points (+/- 0.5, +/- 0.5, +/- 0.5). * ! * ! LATTYP gives the bravais lattice type of the supercell as * ! coded in routine LATGEN (parameter IBRAV) of the * ! lattice vector package (LATTYP may run from 1 to 14). * ! * ! TAUROT,INDEX and WRK are work arrays. * ! * ! NSP is the number of atomic species. * ! NA(NSP) contains the number of atoms per species. * ! * ! NSX and NAX are array dimensioning parameters. * ! * ! PRCHAN contains the unit number for standard output. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! S(3,3,48) contains the INTEGER rotation matrices for all * ! space group operations. * ! * ! SYMOP(3,3,48) contains the rotation matrices of the pure * ! bravais lattice. * ! * ! GTRANS(3,48) contains the nonprimitive translations for * ! all space group operations. * ! * ! NROT contains the number of pure point group rotations (to * ! be stored in S(3,3,1)...S(3,3,NROT)) and * ! NROTK contains the number of space group operations found * ! in routine GETGRP. * ! NOP contains the number of point group operations of the pure * ! bravais lattice without basis. * ! * ! * !*********************************************************************** INTEGER S,SYMOP,SYMGEN(9,3),INV(9),R6Z(9),R3D(9),R2HEX(9) INTEGER R2TRI(9),R2YP(9),R4ZP(9),R4ZBC(9),R4ZFC(9),R2ZP(9) INTEGER R2YBC(9),R2ZBC(9),R2YBAS(9),R2YFC(9),R2ZFC(9),PRCHAN DIMENSION S(3,3,48),GTRANS(3,48),SYMOP(3,3,48),NA(NSP) DIMENSION TAU(NAX,3,NSX),TAUROT(NAX,3,NSX),WRK(NAX,3),INDEX(NAX) SAVE INV,R6Z,R3D,R2HEX,R2TRI,R2YP,R4ZP,R4ZBC,R4ZFC,R2ZP,R2YBC SAVE R2ZBC,R2YBAS,R2YFC,R2ZFC DATA INV /-1,0,0,0,-1,0,0,0,-1/,R3D /0,0,1,1,0,0,0,1,0/ DATA R6Z /1,-1,0,1,0,0,0,0,1/,R2HEX /1,-1,0,0,-1,0,0,0,-1/ DATA R2TRI /-1,0,0,0,0,-1,0,-1,0/,R4ZP /0,-1,0,1,0,0,0,0,1/ DATA R2YP /-1,0,0,0,1,0,0,0,-1/,R4ZBC /0,1,0,0,1,-1,-1,1,0/ DATA R4ZFC /1,1,1,0,0,-1,-1,0,0/,R2ZP /-1,0,0,0,-1,0,0,0,1/ DATA R2YBC /0,-1,1,0,-1,0,1,-1,0/,R2ZBC /0,1,-1,1,0,-1,0,0,-1/ DATA R2YBAS /0,-1,0,-1,0,0,0,0,-1/,R2YFC /0,0,1,-1,-1,-1,1,0,0/ DATA R2ZFC /0,1,0,1,0,0,-1,-1,-1/ ! The pure translation lattice (bravais lattice) has some maximum ! symmetry. Set first up the point group operations for this symmetry. ! Remark: All these groups contain the inversion as a generator ... : CALL SGRCOP(INV,SYMGEN(1,1)) IF (LATTYP==14) THEN ! Triclinic symmetry: CALL SGRGEN(SYMGEN,1,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,14) 14 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' triclinic supercell.'/) ELSE IF (LATTYP==13) THEN ! Monoclinic symmetry (base centered / 'b-a' is rotation axis ...): CALL SGRCOP(R2YBAS,SYMGEN(1,2)) CALL SGRGEN(SYMGEN,2,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,13) 13 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' base centered monoclinic supercell.'/) ELSE IF (LATTYP==12) THEN ! Monoclinic symmetry (primitive cell / 'b-axis' is rotation axis ...): CALL SGRCOP(R2YP,SYMGEN(1,2)) CALL SGRGEN(SYMGEN,2,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,12) 12 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' simple monoclinic supercell.'/) ELSE IF (LATTYP==11) THEN ! Orthorhombic symmetry (base centered ...): CALL SGRCOP(R2ZP,SYMGEN(1,2)) CALL SGRCOP(R2YBAS,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,11) 11 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' base centered orthorhombic supercell.'/) ELSE IF (LATTYP==10) THEN ! Orthorhombic symmetry (face centered ...): CALL SGRCOP(R2ZFC,SYMGEN(1,2)) CALL SGRCOP(R2YFC,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,10) 10 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' face centered orthorhombic supercell.'/) ELSE IF (LATTYP==9) THEN ! Orthorhombic symmetry (body centered ...): CALL SGRCOP(R2ZBC,SYMGEN(1,2)) CALL SGRCOP(R2YBC,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,9) 9 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' body centered orthorhombic supercell.'/) ELSE IF (LATTYP==8) THEN ! Orthorhombic symmetry (primitive cell ...): CALL SGRCOP(R2ZP,SYMGEN(1,2)) CALL SGRCOP(R2YP,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,8) 8 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' simple orthorhombic supercell.'/) ELSE IF (LATTYP==7) THEN ! Trigonal (rhombohedral) symmetry: CALL SGRCOP(R2TRI,SYMGEN(1,2)) CALL SGRCOP(R3D,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,7) 7 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' trigonal (rhombohedral) supercell.'/) ELSE IF (LATTYP==6) THEN ! Tetragonal symmetry (body centered / 'c-axis' is the special axis): CALL SGRCOP(R4ZBC,SYMGEN(1,2)) CALL SGRCOP(R2YBC,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,6) 6 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' body centered tetragonal supercell.'/) ELSE IF (LATTYP==5) THEN ! Tetragonal symmetry (primitive cell / 'c-axis' is the special axis): CALL SGRCOP(R4ZP,SYMGEN(1,2)) CALL SGRCOP(R2YP,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,5) 5 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' simple tetragonal supercell.'/) ELSE IF (LATTYP==4) THEN ! Hexagonal symmetry ('c-axis' is the rotation axis): CALL SGRCOP(R6Z,SYMGEN(1,2)) CALL SGRCOP(R2HEX,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,4) 4 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' hexagonal supercell.'/) ELSE IF (LATTYP==3) THEN ! Cubic symmetry (face centered ...): CALL SGRCOP(R3D,SYMGEN(1,2)) CALL SGRCOP(R4ZFC,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,3) 3 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' face centered cubic supercell.'/) ELSE IF (LATTYP==2) THEN ! Cubic symmetry (body centered ...): CALL SGRCOP(R3D,SYMGEN(1,2)) CALL SGRCOP(R4ZBC,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,2) 2 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' body centered cubic supercell.'/) ELSE IF (LATTYP==1) THEN ! Cubic symmetry (primitive cell ...): CALL SGRCOP(R3D,SYMGEN(1,2)) CALL SGRCOP(R4ZP,SYMGEN(1,3)) CALL SGRGEN(SYMGEN,3,SYMOP,NOP) IF ((PRCHAN>=0).AND.(PRCHAN<=99)) WRITE(PRCHAN,1) 1 FORMAT(/' Routine SETGRP: Setting up the symmetry group for a '/, & & ' simple cubic supercell.'/) END IF ! Now select all symmetry operations which reproduce the lattice: CALL GETGRP(S,GTRANS,NROT,NROTK,SYMOP,NOP,TAU, & & TAUROT,NSP,NSX,NA,NAX,INDEX,WRK,PRCHAN) RETURN END !*********************************************************************** ! * SUBROUTINE GETGRP(S,GTRANS,NROT,NROTK,SYMOP,NOP,TAU, & & TAUROT,NSP,NSX,NA,NAX,INDEX,WRK,PRCHAN) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine GETGRP is a secondary user-interface of this symmetry * ! package. The function of GETGRP is to return all possible space * ! group operators that reproduce a lattice with basis out of a * ! (maximum) pool of point group operations that is compatible with * ! the symmetry of the pure translation lattice without any basis. * ! The setup of the space group operations is done in the following * ! way: We pass through the pool of (possibly allowed) symmetry * ! operations and check each operation whether it can reproduce the * ! lattice with basis - possibly connected with some nontrivial * ! translation. This check will be done by routine CHKSYM (further * ! details will be described there). If the answer is positive the * ! operation will be stored (and counted) - otherwise it will be * ! thrown away. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! SYMOP(3,3,48) contains the rotation matrices of the pure * ! bravais lattice. * ! * ! NOP contains the number of rotation matrices in array SYMOP * ! * ! TAU(NAX,3,NSX) contains the crystal coordinates of all atomic * ! positions within the supercell. Warning: TAU will be * ! modified on output (see comment on routine SETGRP)! * ! * ! TAUROT,INDEX and WRK are work arrays. * ! * ! NSP is the number of atomic species. * ! NA(NSP) contains the number of atoms per species. * ! * ! NSX and NAX are array dimensioning parameters. * ! * ! PRCHAN contains the unit number for standard output. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! S(3,3,48) contains the INTEGER rotation matrices for all * ! space group operations. * ! * ! GTRANS(3,48) contains the nonprimitive translations for * ! all space group operations. * ! * ! NROT contains the number of pure point group rotations (to * ! be stored in S(3,3,1)...S(3,3,NROT)) and * ! NROTK contains the number of space group operations found * ! in routine GETGRP. * ! * ! * !*********************************************************************** LOGICAL TS INTEGER S,SYMOP,ZERO(9),HELP(3,3,48),PRCHAN DIMENSION S(3,3,48),GTRANS(3,48),SYMOP(3,3,48),TEMP(3,48),NA(NSP) DIMENSION TAU(NAX,3,NSX),TAUROT(NAX,3,NSX),WRK(NAX,3),INDEX(NAX) SAVE TINY,ZERO DATA TINY / MY_TINY /,ZERO /0,0,0,0,0,0,0,0,0/ NROT=0 NROTK=0 ! Check the operations of a 'pool' of all possible operations ... : DO 1 I=1,NOP CALL CHKSYM(SYMOP(1,1,I),GTRANS(1,I),TS,TAU,TAUROT, & & NSP,NSX,NA,NAX,INDEX,WRK) IF (TS) THEN ! Hurray! We found some symmetry operation in subroutine CHKSYM ... : IF ((ABS(GTRANS(1,I))0) THEN ! If there are operations with nontrivial translations then store them ! into the array S and GTRANS (placing them after the last pure point ! group operations S(I,J,NROT) ): DO 2 I=1,NROTK CALL SGRCOP(HELP(1,1,I),S(1,1,I+NROT)) GTRANS(1,I+NROT)=TEMP(1,I) GTRANS(2,I+NROT)=TEMP(2,I) GTRANS(3,I+NROT)=TEMP(3,I) 2 CONTINUE END IF ! Total number of space group operations: NROTK=NROTK+NROT IF (NROTK<48) THEN ! Fill the rest of arrays S and GTRANS with zeros ... : DO 3 I=NROTK+1,48 CALL SGRCOP(ZERO,S(1,1,I)) GTRANS(1,I)=0._q GTRANS(2,I)=0._q GTRANS(3,I)=0._q 3 CONTINUE END IF IF ((PRCHAN>=0).AND.(PRCHAN<=99)) & & WRITE(PRCHAN,100) NROTK,NROT,NOP 100 FORMAT(/' Subroutine GETGRP returns: Found ',I2, & & ' space group operations'/,' (whereof ',I2, & & ' operations were pure point group operations)'/, & & ' out of a pool of ',I2,' trial point group operations.'/) RETURN END !*********************************************************************** ! * SUBROUTINE CHKSYM(S,GTRANS,TS,TAU,TAUROT,NSP,NSX,NA,NAX,INDEX,WRK) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine CHKSYM checks whether a point group symmetry element is * ! a valid symmetry operation on a supercell (possibly connected with * ! some nontrivial translation) or not. The procedure is in principle * ! very easy: The point group transformation will be imposed on all * ! atomic positions (for all species) and then the rotated supercell * ! will be compared with the original cell. There are only three * ! possibilities: All positions will be reproduced (pure point group * ! operation), all positions will be reproduced after an additional * ! unique translation of all atoms (general space group operation) * ! or it will not be possible to reproduce the supercell in any way * ! (no allowed operation). The translation can be guessed quite easy * ! by looking on the first atom: If we try all difference vectors * ! between the position of the first atom of the unrotated cell and * ! the positions of all atoms of the rotated cell then one of these * ! difference vectors must be the translation vector (if it exists!). * ! This vector must apply to all atoms. So we translate all atoms in * ! the rotated cell by minus this vector - this should reproduce the * ! original positions of the unrotated cell. The central point is now * ! only how to compare the rotated and the unrotated supercell. The * ! algorithm used bases on the following idea: If we sort all atomic * ! positions within a cell as it is described in routine LATORD of * ! the lattice vector package (sorting first the x-coordinates, then * ! sorting the y-coordinates of all atoms with equal x-coordinates * ! and finally sorting all z-cordinates for all atomic positions with * ! equal x- and y-coordinates - or using other ordering sequences * ! instead of xyz-ordering ...) then we get some unique sequence of * ! ordered positions. Because the same set of input data will mean * ! identical output we can simply perform a one-by-one comparison of * ! the unrotated (ordered) and the rotated (ordered and translated) * ! coordinates. If we have found some group operation the comparison * ! will give the value true ('all positions are equal'). In all other * ! cases the result will be false ('some positions are not equal'! * ! There is just one point one has to care about: Supercells which * ! are non-primitive. Then there exist also 'non-primitive primitive' * ! translations (trivial translations of the generating cell ...)! * ! To avoid trouble with non-uniquely defined translations in this * ! case we just take the vectors with smallest coordinates (from all * ! the vectors which will be detected to be an allowed translation). * ! * ! It should be trivial but will be remarked finally: The procedure * ! must be performed on every atomic species seperately and the * ! result (allowed operation or not / nontrivial translation) must be * ! the same for all atomic species - otherwise the symmetry operation * ! is not a valid operation ... . * ! * ! Final very important notice: Comparison of atomic positions will * ! be done by looking at differences of coordinates (which must be * ! smaller than some tolerance). The allowed tolerances for the * ! coordinates are very small (1E-9 absolute)!!! Therefore all input * ! coordinates should be given with an accuray of at least +/-1E-10 * ! because otherwise this routine could work incorrectly (could fail * ! to detect space group operations ...)! * ! * ! * ! Input parameters: * ! ----------------- * ! * ! S(3,3) is the point group symmetry element to be checked for * ! validity (INTEGER rotation matrix). * ! * ! TAU(NAX,3,NSX) contains the crystal coordinates of all atomic * ! positions within the supercell. Warning: TAU will be * ! modified on output (see comment on routine SETGRP)! * ! * ! TAUROT,INDEX and WRK are work arrays. * ! * ! NSP is the number of atomic species. * ! NA(NSP) contains the number of atoms per species. * ! * ! NSX and NAX are array dimensioning parameters. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! GTRANS(3) contains the nonprimitive translations which is * ! connected with the rotation matrix S(3,3). * ! * ! TS is a logical value which returns whether S(3,3) is a valid * ! symmetry operation on the supercell or not. * ! * ! * !*********************************************************************** LOGICAL TS,TNTRAN,TCONT INTEGER S DIMENSION S(3,3),GTRANS(3),NA(NSP),TAU(NAX,3,NSX),TAUDIF(3) DIMENSION TAUROT(NAX,3,NSX),WRK(NAX,3),INDEX(NAX),TAUSAV(3),TRA(3) SAVE TINY DATA TINY / MY_TINY / ! Unfortunately one has always two possibilities for the book-keeping ! of atomic coordinates, forces etc.: One can distinguish between the ! species and the atoms within each species (two dimensional arrays ! "NAX,NSX") or one can only count the atoms (one dimensional arrays ! "NSAX") where the species are distinguished be the map NA(NSP) in ! the sense (and this should be the convention used - otherwise this ! symmetry package will not work correctly!!) that the atoms of species ! "1" are stored in elements 1...NA(1), atoms of species "2" in elements ! NA(1)+1...NA(1)+NA(2) and so on and atoms of species "NSP" are stored ! in elements NA(1)+NA(2)+...+NA(NSP-1)+1 ... NA(1)+NA(2)+...+NA(NSP). ! The package works for both cases - in the second case one must set ! the dimensioning parameter NSX=1 (only important for NSP>1) ... . TCONT=((NSX==1).AND.(NSP>1)) ! Now start ... : ISMIN=1 TRA(1)=2._q TRA(2)=2._q TRA(3)=2._q TS=.FALSE. ISTART=1 IMINST=1 ! For all atomic species ... : DO 3 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) ! Search the species with smallest number of atoms ... : IF (NA(IS) [-0.5,0.5) ): TAU(IA,1,ISP)=TAU(IA,1,ISP)-NINT(TAU(IA,1,ISP)) TAU(IA,2,ISP)=TAU(IA,2,ISP)-NINT(TAU(IA,2,ISP)) TAU(IA,3,ISP)=TAU(IA,3,ISP)-NINT(TAU(IA,3,ISP)) TAU(IA,1,ISP)=MOD(TAU(IA,1,ISP)+6.5_q,1._q)-0.5_q TAU(IA,2,ISP)=MOD(TAU(IA,2,ISP)+6.5_q,1._q)-0.5_q TAU(IA,3,ISP)=MOD(TAU(IA,3,ISP)+6.5_q,1._q)-0.5_q ! It could happen that (due to roundoff errors!) the surface coordinates ! are not uniquely -0.5 but probably +0.49999999999... . Ensure that all ! these coordinates are transformed to -0.5! IF ((ABS(TAU(IA,1,ISP)+0.5_q)) [-0.5,0.5) ): TAUROT(IA,1,ISP)=MOD(TAUROT(IA,1,ISP)+6.5_q,1._q)-0.5_q TAUROT(IA,2,ISP)=MOD(TAUROT(IA,2,ISP)+6.5_q,1._q)-0.5_q TAUROT(IA,3,ISP)=MOD(TAUROT(IA,3,ISP)+6.5_q,1._q)-0.5_q ! It could happen that (due to roundoff errors!) the surface coordinates ! are not uniquely -0.5 but probably +0.49999999999... . Ensure that all ! these coordinates are transformed to -0.5! IF ((ABS(TAUROT(IA,1,ISP)+0.5_q))(TRA(1)+TINY)).OR. & & (GTRANS(2)>(TRA(2)+TINY)).OR. & & (GTRANS(3)>(TRA(3)+TINY))) GOTO 10 ! Translate the rotated coordinates by GTRANS ... ISTART=1 DO 5 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) DO 4 IA=ISTART,NA(IS)+ISTART-1 TAUROT(IA,1,ISP)=TAUROT(IA,1,ISP)+GTRANS(1) TAUROT(IA,2,ISP)=TAUROT(IA,2,ISP)+GTRANS(2) TAUROT(IA,3,ISP)=TAUROT(IA,3,ISP)+GTRANS(3) ! ... impose the periodic boundary condition on these coordinates ... TAUROT(IA,1,ISP)=MOD(TAUROT(IA,1,ISP)+6.5_q,1._q)-0.5_q TAUROT(IA,2,ISP)=MOD(TAUROT(IA,2,ISP)+6.5_q,1._q)-0.5_q TAUROT(IA,3,ISP)=MOD(TAUROT(IA,3,ISP)+6.5_q,1._q)-0.5_q ! ... ensuring that all surface coordinates are put to -0.5 ... IF ((ABS(TAUROT(IA,1,ISP)+0.5_q))1)) ISTART=ISTART+NA(IS-1) DO 6 IA=ISTART,NA(IS)+ISTART-1 ! Take the difference of the rotated and the unrotated coordinates after ! we translated them - the result should be (0.,0.,0.)! TAUDIF(1)=TAU(IA,1,ISP)-TAUROT(IA,1,ISP) TAUDIF(2)=TAU(IA,2,ISP)-TAUROT(IA,2,ISP) TAUDIF(3)=TAU(IA,3,ISP)-TAUROT(IA,3,ISP) ! TAUDIF may only differ from (0.,0.,0.) by some trivial translation: TAUDIF(1)=MOD((TAUDIF(1)+6._q),1._q) TAUDIF(2)=MOD((TAUDIF(2)+6._q),1._q) TAUDIF(3)=MOD((TAUDIF(3)+6._q),1._q) ! For reasons of safety: IF ((ABS(TAUDIF(1)-1._q))<(TINY*0.5_q)) TAUDIF(1)=0._q IF ((ABS(TAUDIF(2)-1._q))<(TINY*0.5_q)) TAUDIF(2)=0._q IF ((ABS(TAUDIF(3)-1._q))<(TINY*0.5_q)) TAUDIF(3)=0._q ! If TAUDIF is not equal (0.,0.,0.) TNTRAN=TNTRAN.AND.FALSE and so it ! remains FALSE forever ... (if only one position is not reproduced!). ! Only if all TAUDIFs are zero TNTRAN (starting with the value TRUE) ! will remain TRUE (that means we found an allowed symmetry operation). TNTRAN=TNTRAN.AND.((ABS(TAUDIF(1))1)) ISTART=ISTART+NA(IS-1) DO 8 IA=ISTART,NA(IS)+ISTART-1 TAUROT(IA,1,ISP)=TAUROT(IA,1,ISP)-GTRANS(1) TAUROT(IA,2,ISP)=TAUROT(IA,2,ISP)-GTRANS(2) TAUROT(IA,3,ISP)=TAUROT(IA,3,ISP)-GTRANS(3) 8 CONTINUE 9 CONTINUE 10 CONTINUE IF (TS) THEN GTRANS(1)=TRA(1) GTRANS(2)=TRA(2) GTRANS(3)=TRA(3) END IF ! Uff ... . Bye bye! RETURN END !*********************************************************************** ! * SUBROUTINE DYNSYM(V,MAP,S,T,NROTK,NROT,PT,NP, & & NSX,NSP,NAX,NA,NSAX,VR,W,IND,P) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine DYNSYM tests whether an array of vectors (e.g. velocities) * ! located at different atomic positions is compatible with a given * ! symmetry operation. The test is done by symmetrizing the vectors * ! as in routine FSYM (see above) and comparing the result to the * ! non-symmetrized vectors. If the test will be negative the symmetry * ! operation will be discarded from the set of given operations. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! V(NAX,3,NSX) contains the crystal coordinates of the vectors * ! located at the corresponding atomic position. * ! * ! MAP(NAX,NSX,48) contains a table connecting the rotated and * ! the unrotated positions (stored is the index of * ! the rotated position = first index in array TAU). * ! MAP can be generated by routine CLASS (see below). * ! S(3,3,48) contains the INTEGER rotation matrices. * ! T(3,48) contains the nonprimitive translations. * ! NROTK contains the number of given space group operations. * ! NROT contains the number of pure point group operations. * ! PT contains all "primitive" translations associated with the * ! generating unit cell of a nonprimitive supercell * ! NP is the number of "primitive" translations in PT * ! * ! NSP is the number of atomic species. * ! NA(NSP) contains the number of atoms per species. * ! * ! NSX, NAX and NSAX are array dimensioning parameters. * ! * ! VR, W and IND are work arrays. * ! * ! P contains the unit number for standard output. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! Arrays S and T contain the compatible operations on output * ! (input is overwritten!), NROTK and NROT the numbers of elements * ! in S and T being space group / point group operations. * ! * ! * !*********************************************************************** INTEGER S,GRP,P,OLDNR,OLDNRK LOGICAL TEST,TCONT DIMENSION V(NAX,3,NSX),MAP(NAX,NSX,48,NP),S(3,3,48),NA(NSP) DIMENSION VR(NAX,3),GRP(3,3,48),T(3,48),TCOP(3,48),IM(48) DIMENSION PT(NSAX+2,3),IND(NSAX),W(NAX,3) SAVE TINY DATA TINY / MY_TINY / ! Unfortunately one has always two possibilities for the book-keeping ! of atomic coordinates, forces etc.: One can distinguish between the ! species and the atoms within each species (two dimensional arrays ! "NAX,NSX") or one can only count the atoms (one dimensional arrays ! "NSAX") where the species are distinguished be the map NA(NSP) in ! the sense (and this should be the convention used - otherwise this ! symmetry package will not work correctly!!) that the atoms of species ! "1" are stored in elements 1...NA(1), atoms of species "2" in elements ! NA(1)+1...NA(1)+NA(2) and so on and atoms of species "NSP" are stored ! in elements NA(1)+NA(2)+...+NA(NSP-1)+1 ... NA(1)+NA(2)+...+NA(NSP). ! The package works for both cases - in the second case one must set ! the dimensioning parameter NSX=1 (only important for NSP>1) ... . TCONT=((NSX==1).AND.(NSP>1)) ! Now start ... : OLDNR=NROT OLDNRK=NROTK DO 1 IROT=1,48 1 IM(IROT)=IROT ISTART=1 ! For all atomic species ... : DO 7 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) ! First copy the rotation matrices and translation vectors: DO 2 IROT=1,144 GRP(1,IROT,1)=S(1,IROT,1) GRP(2,IROT,1)=S(2,IROT,1) GRP(3,IROT,1)=S(3,IROT,1) TCOP(IROT,1)=T(IROT,1) S(1,IROT,1)=0 S(2,IROT,1)=0 S(3,IROT,1)=0 T(IROT,1)=0 2 CONTINUE ! Copy the input vectors: DO 3 IA=ISTART,NA(IS)+ISTART-1 W(IA,1)=V(IA,1,ISP) W(IA,2)=V(IA,2,ISP) W(IA,3)=V(IA,3,ISP) 3 CONTINUE NR=0 ! Apply all space group operations on the vectors: DO 6 IROT=1,NROTK ! Calculate the rotated vectors for all atoms: DO 4 IA=ISTART,NA(IS)+ISTART-1 VR(IA,1)=GRP(1,1,IROT)*W(IA,1)+ & & GRP(2,1,IROT)*W(IA,2)+ & & GRP(3,1,IROT)*W(IA,3) VR(IA,2)=GRP(1,2,IROT)*W(IA,1)+ & & GRP(2,2,IROT)*W(IA,2)+ & & GRP(3,2,IROT)*W(IA,3) VR(IA,3)=GRP(1,3,IROT)*W(IA,1)+ & & GRP(2,3,IROT)*W(IA,2)+ & & GRP(3,3,IROT)*W(IA,3) 4 CONTINUE ! Compare them with the unrotated vectors: TEST=.TRUE. DO 5 IA=ISTART,NA(IS)+ISTART-1 ! TEST will remain true if all rotated vectors are equal to the original ! vectors; otherwise TEST will become false! JA=MAP(IA,ISP,IM(IROT),1) TEST=TEST.AND. & & (ABS(V(IA,1,ISP)-VR(JA,1))1)) ISTART=ISTART+NA(IS-1) ! First copy all translation vectors: DO 9 ITRANS=1,NP W(ITRANS,1)=PT(ITRANS,1) W(ITRANS,2)=PT(ITRANS,2) W(ITRANS,3)=PT(ITRANS,3) PT(ITRANS,1)=0._q PT(ITRANS,2)=0._q PT(ITRANS,3)=0._q 9 CONTINUE NPRI=0 TEST=.TRUE. ! Apply all primitive translations: DO 11 ITRANS=1,NP ! Compare vectors with vectors at translated positions: DO 10 IA=ISTART,NA(IS)+ISTART-1 ! TEST will remain true if all tranlated vectors are equal to the ! original vectors; otherwise TEST will become false! JA=MAP(IA,ISP,1,IND(ITRANS)) TEST=TEST.AND. & & (ABS(V(IA,1,ISP)-V(JA,1,ISP))=0).AND.(P<=99)) WRITE(P,50) NROTK,NROT,OLDNRK,OLDNR,NP 50 FORMAT(/' Subroutine DYNSYM returns: Found ',I2, & & ' space group operations'/,' (whereof ',I2, & & ' operations were pure point group operations)'/, & & ' out of a pool of ',I2,' trial space group operations'/, & & ' (whereof ',I2, & & ' operations were pure point group operations)'/, & & ' and found also ',I5,' ''primitive'' translations'/) RETURN END !*********************************************************************** ! * SUBROUTINE MAGSYM(ATOMOM,MAGROT,MAP,S,T,NROTK,NROT,PT,NP, & & NSX,NSP,NAX,NA,NSAX,MAGR,W,IND,P) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine MAGSYM tests whether a given pattern of magnetic moments * ! located at different atomic positions is compatible with a given * ! symmetry operation. The test is done by checking atomic magnetic * ! moments at symmetry-related atoms -- they can be identical or may * ! differ by a factor (-1). If the test will be negative the symmetry * ! operation will be discarded from the set of given operations. If * ! the test was positive the connection factor (1 or -1) is stored. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! ATOMOM(NAX,NSX) contains the atomic magnetic moments located * ! at the corresponding atomic position. * ! * ! MAP(NAX,NSX,48) contains a table connecting the rotated and * ! the unrotated positions (stored is the index of * ! the rotated position = first index in array TAU). * ! MAP can be generated by routine CLASS (see below). * ! S(3,3,48) contains the INTEGER rotation matrices. * ! T(3,48) contains the nonprimitive translations. * ! NROTK contains the number of given space group operations. * ! NROT contains the number of pure point group operations. * ! PT contains all "primitive" translations associated with the * ! generating unit cell of a nonprimitive supercell * ! NP is the number of "primitive" translations in PT * ! * ! NSP is the number of atomic species. * ! NA(NSP) contains the number of atoms per species. * ! * ! NSX, NAX and NSAX are array dimensioning parameters. * ! * ! MAGR, W and IND are work arrays. * ! * ! P contains the unit number for standard output. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! Array MAGMOM(48,NP) contains the connection factors between a * ! magnetic moment pattern and the pattern at the rotated atomic * ! positions. In the ISING case we have two possibilities: 1 or -1 * ! * ! Arrays S and T contain the compatible operations on output * ! (input is overwritten!), NROTK and NROT the numbers of elements * ! in S and T being space group / point group operations. * ! * ! * !*********************************************************************** REAL(q) MAGROT,MAGR,MR INTEGER S,GRP,P,OLDNR,OLDNRK LOGICAL TEST,TCONT DIMENSION MAGROT(48,NSAX),MAP(NAX,NSX,48,NP),S(3,3,48),NA(NSP) DIMENSION MAGR(NAX,3),GRP(3,3,48),T(3,48),TCOP(3,48),IM(48) DIMENSION PT(NSAX+2,3),IND(NSAX),W(NAX,3),ATOMOM(NAX,NSX) SAVE TINY DATA TINY / MY_TINY / ! Unfortunately one has always two possibilities for the book-keeping ! of atomic coordinates, forces etc.: One can distinguish between the ! species and the atoms within each species (two dimensional arrays ! "NAX,NSX") or one can only count the atoms (one dimensional arrays ! "NSAX") where the species are distinguished be the map NA(NSP) in ! the sense (and this should be the convention used - otherwise this ! symmetry package will not work correctly!!) that the atoms of species ! "1" are stored in elements 1...NA(1), atoms of species "2" in elements ! NA(1)+1...NA(1)+NA(2) and so on and atoms of species "NSP" are stored ! in elements NA(1)+NA(2)+...+NA(NSP-1)+1 ... NA(1)+NA(2)+...+NA(NSP). ! The package works for both cases - in the second case one must set ! the dimensioning parameter NSX=1 (only important for NSP>1) ... . TCONT=((NSX==1).AND.(NSP>1)) ! Now start ... : OLDNR=NROT OLDNRK=NROTK DO 1 IROT=1,48 1 IM(IROT)=IROT ! Initialize array MAGMOM: DO 2 IROT=1,48*NSAX MAGROT(IROT,1)=0._q 2 CONTINUE ! The first task will be to search the first atom with a non-vanishing moment TEST=.TRUE. IAMAG=0 ISMAG=0 ISTART=1 DO 4 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) DO 3 IA=ISTART,NA(IS)+ISTART-1 TEST=TEST.AND.(ABS(ATOMOM(IA,ISP)) no symmetry breaking! 5 IF (TEST) THEN DO 7 ITRANS=1,NP DO 6 IROT=1,NROTK MAGROT(IROT,ITRANS)=1._q 6 CONTINUE 7 CONTINUE GOTO 19 ENDIF ISTART=1 ! For all atomic species ... : DO 12 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) ! First copy the rotation matrices and translation vectors: DO 8 IROT=1,144 GRP(1,IROT,1)=S(1,IROT,1) GRP(2,IROT,1)=S(2,IROT,1) GRP(3,IROT,1)=S(3,IROT,1) TCOP(IROT,1)=T(IROT,1) S(1,IROT,1)=0 S(2,IROT,1)=0 S(3,IROT,1)=0 T(IROT,1)=0 8 CONTINUE NR=0 ! Apply all space group operations to the atomic coordinates: DO 11 IROT=1,NROTK ! Rotate the first atom of the first species having a non-zero moment: JA=MAP(IAMAG,ISMAG,IM(IROT),1) TEST =(ABS(ABS(ATOMOM(IAMAG,ISMAG))-ABS(ATOMOM(JA,ISMAG)))1)) ISTART=ISTART+NA(IS-1) ! First copy all translation vectors: DO 14 ITRANS=1,NP W(ITRANS,1)=PT(ITRANS,1) W(ITRANS,2)=PT(ITRANS,2) W(ITRANS,3)=PT(ITRANS,3) PT(ITRANS,1)=0._q PT(ITRANS,2)=0._q PT(ITRANS,3)=0._q 14 CONTINUE NPRI=0 TEST=.TRUE. ! Apply all primitive translations: DO 17 ITRANS=1,NP ! Translate the first atom of the first species having a non-zero moment: JA=MAP(IAMAG,ISMAG,1,IND(ITRANS)) TEST =(ABS(ABS(ATOMOM(IAMAG,ISMAG))-ABS(ATOMOM(JA,ISMAG)))=0).AND.(P<=99)) WRITE(P,50) NROTK,NROT,OLDNRK,OLDNR,NP ! And now the most interesting moment: I claim that combined application of ! a rotation IROT and a translation ITRANS should result in a multiplication ! of the corresponding transformations for the magnetic moment pattern ... : DO 21 IROT=2,NROTK DO 20 ITRANS=2,NP MAGROT(IROT,ITRANS)=MAGROT(IROT,1)*MAGROT(1,ITRANS) ! is this true? 20 CONTINUE 21 CONTINUE 50 FORMAT(/' Subroutine MAGSYM returns: Found ',I2, & & ' space group operations'/,' (whereof ',I2, & & ' operations were pure point group operations)'/, & & ' out of a pool of ',I2,' trial space group operations'/, & & ' (whereof ',I2, & & ' operations were pure point group operations)'/, & & ' and found also ',I5,' ''primitive'' translations'/) RETURN END !*********************************************************************** ! * SUBROUTINE FSYM(F,MAP,S,NR,NP,NSX,NSP,NAX,NA,FR,W, & & A1,A2,A3,B1,B2,B3) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine FSYM symmetrizes an array of vectors (forces, velocities) * ! located at different atomic positions. The procedure is a quite * ! simple: Rotate the input vector and add it to the vector at the * ! rotated atomic position (for all space group operations ...). * ! Finally divide the result by the number of symmetry operations! * ! * ! * ! Input parameters: * ! ----------------- * ! * ! F(3,NAX,NSX) contains the cartesian coordinates of the vector * ! located at the corresponding atomic position. * ! * ! MAP(NAX,NSX,48,NP) contains a table connecting the rotated and * ! the unrotated positions (stored is the index of * ! the rotated position = first index in array TAU). * ! MAP can be generated by routine CLASS (see below). * ! S(3,3,48) contains the INTEGER rotation matrices. * ! NR contains the number of given space group operations. * ! NP contains the number of "primitive" translations in the cell. * ! * ! NSP is the number of atomic species. * ! NA(NSP) contains the number of atoms per species. * ! * ! NSX and NAX are array dimensioning parameters. * ! * ! FR and W are work arrays. * ! * ! A1, A2 and A3 contain the direct (real space) lattice vectors. * ! B1, B2 and B3 contain the reciprocal lattice vectors. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! Array F contains the symmetrized vectors on output (input is * ! overwritten!). * ! * ! * !*********************************************************************** LOGICAL TCONT INTEGER S DIMENSION F(3,NAX,NSX),MAP(NAX,NSX,48,NP),S(3,3,48),NA(NSP) DIMENSION FR(NAX,3),W(NAX,3),A1(3),A2(3),A3(3),B1(3),B2(3),B3(3) ! Unfortunately one has always two possibilities for the book-keeping ! of atomic coordinates, forces etc.: One can distinguish between the ! species and the atoms within each species (two dimensional arrays ! "NAX,NSX") or one can only count the atoms (one dimensional arrays ! "NSAX") where the species are distinguished be the map NA(NSP) in ! the sense (and this should be the convention used - otherwise this ! symmetry package will not work correctly!!) that the atoms of species ! "1" are stored in elements 1...NA(1), atoms of species "2" in elements ! NA(1)+1...NA(1)+NA(2) and so on and atoms of species "NSP" are stored ! in elements NA(1)+NA(2)+...+NA(NSP-1)+1 ... NA(1)+NA(2)+...+NA(NSP). ! The package works for both cases - in the second case one must set ! the dimensioning parameter NSX=1 (only important for NSP>1) ... . TCONT=((NSX==1).AND.(NSP>1)) ! Now start ... : ISTART=1 ! For all atomic species ... : DO 7 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) ! First transform the cartesian coordiates to crystal coordinates ... : DO 1 IA=ISTART,NA(IS)+ISTART-1 W(IA,1)=B1(1)*F(1,IA,ISP)+B1(2)*F(2,IA,ISP)+ & & B1(3)*F(3,IA,ISP) W(IA,2)=B2(1)*F(1,IA,ISP)+B2(2)*F(2,IA,ISP)+ & & B2(3)*F(3,IA,ISP) W(IA,3)=B3(1)*F(1,IA,ISP)+B3(2)*F(2,IA,ISP)+ & & B3(3)*F(3,IA,ISP) 1 CONTINUE ! Reset array F and divide the input vectors by NR*NP: DO 2 IA=ISTART,NA(IS)+ISTART-1 W(IA,1)=W(IA,1)/FLOAT(NR*NP) W(IA,2)=W(IA,2)/FLOAT(NR*NP) W(IA,3)=W(IA,3)/FLOAT(NR*NP) F(1,IA,ISP)=0._q F(2,IA,ISP)=0._q F(3,IA,ISP)=0._q 2 CONTINUE ! Apply all space group operations on the vectors: DO 6 IROT=1,NR ! Calculate the rotated vectors for all atoms: DO 3 IA=ISTART,NA(IS)+ISTART-1 FR(IA,1)=S(1,1,IROT)*W(IA,1)+ & & S(2,1,IROT)*W(IA,2)+ & & S(3,1,IROT)*W(IA,3) FR(IA,2)=S(1,2,IROT)*W(IA,1)+ & & S(2,2,IROT)*W(IA,2)+ & & S(3,2,IROT)*W(IA,3) FR(IA,3)=S(1,3,IROT)*W(IA,1)+ & & S(2,3,IROT)*W(IA,2)+ & & S(3,3,IROT)*W(IA,3) 3 CONTINUE ! Apply all translations to all atoms: DO 5 ITRANS=1,NP ! Now add all rotated vectors to the corresponding atomic vector: DO 4 IA=ISTART,NA(IS)+ISTART-1 F(1,IA,ISP)=F(1,IA,ISP)+ & & FR(MAP(IA,ISP,IROT,ITRANS),1)*A1(1)+ & & FR(MAP(IA,ISP,IROT,ITRANS),2)*A2(1)+ & & FR(MAP(IA,ISP,IROT,ITRANS),3)*A3(1) F(2,IA,ISP)=F(2,IA,ISP)+ & & FR(MAP(IA,ISP,IROT,ITRANS),1)*A1(2)+ & & FR(MAP(IA,ISP,IROT,ITRANS),2)*A2(2)+ & & FR(MAP(IA,ISP,IROT,ITRANS),3)*A3(2) F(3,IA,ISP)=F(3,IA,ISP)+ & & FR(MAP(IA,ISP,IROT,ITRANS),1)*A1(3)+ & & FR(MAP(IA,ISP,IROT,ITRANS),2)*A2(3)+ & & FR(MAP(IA,ISP,IROT,ITRANS),3)*A3(3) 4 CONTINUE 5 CONTINUE 6 CONTINUE 7 CONTINUE ! Ready ... ! RETURN END !*********************************************************************** ! * SUBROUTINE PSYM(LPBC,X,MAP,S,GT,NR,PT,NP,NSX,NSP,NAX,NA,XR,W, & & A1,A2,A3,B1,B2,B3) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine PSYM symmetrizes atomic coordinates. The procedure is * ! almost the same as in routine FSYM. The only differences are: * ! We may not only rotate the vectors X (playing a similar role * ! like F in routine FSYM). We must also apply all translations and * ! periodic boundary conditions (= appropriate triviav translations. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! LPBC determines whether (if set .TRUE.) to keep the periodic * ! boundary conditions imposed on the coordinates or * ! not (by trivial back-translations if set .FALSE.). * ! * ! X(3,NAX,NSX) contains the cartesian coordinates of the atoms. * ! * ! MAP(NAX,NSX,48,NP) contains a table connecting the rotated and * ! the unrotated positions (stored is the index of * ! the rotated position = first index in array TAU). * ! MAP can be generated by routine CLASS (see below). * ! S(3,3,48) contains the INTEGER rotation matrices. * ! GT(3,48) contains the corresponding non-trivial translations. * ! NR contains the number of given space group operations. * ! PT(NAX+2,3) contains all "primitive" translations in the cell. * ! NP contains the number of "primitive" translations in the cell. * ! * ! NSP is the number of atomic species. * ! NA(NSP) contains the number of atoms per species. * ! * ! NSX and NAX are array dimensioning parameters. * ! * ! XR and W are work arrays. * ! * ! A1, A2 and A3 contain the direct (real space) lattice vectors. * ! B1, B2 and B3 contain the reciprocal lattice vectors. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! Array X contains the symmetrized vectors on output (input is * ! overwritten!). * ! * ! * !*********************************************************************** LOGICAL LPBC,TCONT INTEGER S DIMENSION X(3,NAX,NSX),MAP(NAX,NSX,48,NP),S(3,3,48),NA(NSP) DIMENSION XR(NAX,3),W(NAX,3),A1(3),A2(3),A3(3),B1(3),B2(3),B3(3) DIMENSION GT(3,48),PT(NAX+2,3) SAVE TINY DATA TINY / MY_TINY / ! Unfortunately one has always two possibilities for the book-keeping ! of atomic coordinates, forces etc.: One can distinguish between the ! species and the atoms within each species (two dimensional arrays ! "NAX,NSX") or one can only count the atoms (one dimensional arrays ! "NSAX") where the species are distinguished be the map NA(NSP) in ! the sense (and this should be the convention used - otherwise this ! symmetry package will not work correctly!!) that the atoms of species ! "1" are stored in elements 1...NA(1), atoms of species "2" in elements ! NA(1)+1...NA(1)+NA(2) and so on and atoms of species "NSP" are stored ! in elements NA(1)+NA(2)+...+NA(NSP-1)+1 ... NA(1)+NA(2)+...+NA(NSP). ! The package works for both cases - in the second case one must set ! the dimensioning parameter NSX=1 (only important for NSP>1) ... . TCONT=((NSX==1).AND.(NSP>1)) ! Now start ... : ISTART=1 ! For all atomic species ... : DO 7 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) ! First transform the cartesian coordiates to crystal coordinates ... : DO 1 IA=ISTART,NA(IS)+ISTART-1 W(IA,1)=B1(1)*X(1,IA,ISP)+B1(2)*X(2,IA,ISP)+ & & B1(3)*X(3,IA,ISP) W(IA,2)=B2(1)*X(1,IA,ISP)+B2(2)*X(2,IA,ISP)+ & & B2(3)*X(3,IA,ISP) W(IA,3)=B3(1)*X(1,IA,ISP)+B3(2)*X(2,IA,ISP)+ & & B3(3)*X(3,IA,ISP) 1 CONTINUE ! Preset array X: DO 2 IA=ISTART,NA(IS)+ISTART-1 IF (.NOT.LPBC) THEN ! Copy W to X X(1,IA,ISP)=W(IA,1) X(2,IA,ISP)=W(IA,2) X(3,IA,ISP)=W(IA,3) ! Apply peridodic boundary conditions to X, i.e. shift back to [-0.5,0.5): X(1,IA,ISP)=MOD(X(1,IA,ISP)+60.5_q,1._q)-0.5_q X(2,IA,ISP)=MOD(X(2,IA,ISP)+60.5_q,1._q)-0.5_q X(3,IA,ISP)=MOD(X(3,IA,ISP)+60.5_q,1._q)-0.5_q ! For reasons of safety (atoms very close to the border of the unit cell): IF (ABS(X(1,IA,ISP)-0.5_q)1) ... . TCONT=((NSX==1).AND.(NSP>1)) ! Now start ... : ISTART=1 ! For all atomic species ... : DO 7 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) DO 1 IA=ISTART,NA(IS)+ISTART-1 ! First impose periodic boundary conditions on TAU (--> [-0.5,0.5) ): TAU(IA,1,ISP)=TAU(IA,1,ISP)-INT(TAU(IA,1,ISP)) TAU(IA,2,ISP)=TAU(IA,2,ISP)-INT(TAU(IA,2,ISP)) TAU(IA,3,ISP)=TAU(IA,3,ISP)-INT(TAU(IA,3,ISP)) TAU(IA,1,ISP)=MOD(TAU(IA,1,ISP)+6.5_q,1._q)-0.5_q TAU(IA,2,ISP)=MOD(TAU(IA,2,ISP)+6.5_q,1._q)-0.5_q TAU(IA,3,ISP)=MOD(TAU(IA,3,ISP)+6.5_q,1._q)-0.5_q ! It could happen that (due to roundoff errors!) the surface coordinates ! are not uniquely -0.5 but probably +0.49999999999... . Ensure that all ! these coordinates are transformed to -0.5! IF ((ABS(TAU(IA,1,ISP)+0.5_q))1) ... . TCONT=((NSX==1).AND.(NSP>1)) ! Now start ... : ISTART=1 ! For all atomic species ... : DO 8 IS=1,NSP ISP=IS IF (TCONT) ISP=1 IF (TCONT.AND.(IS>1)) ISTART=ISTART+NA(IS-1) NCL(IS)=0 DO 1 IA=ISTART,NA(IS)+ISTART-1 ! First impose periodic boundary conditions on TAU (--> [-0.5,0.5) ): TAU(IA,1,ISP)=TAU(IA,1,ISP)-INT(TAU(IA,1,ISP)) TAU(IA,2,ISP)=TAU(IA,2,ISP)-INT(TAU(IA,2,ISP)) TAU(IA,3,ISP)=TAU(IA,3,ISP)-INT(TAU(IA,3,ISP)) TAU(IA,1,ISP)=MOD(TAU(IA,1,ISP)+6.5_q,1._q)-0.5_q TAU(IA,2,ISP)=MOD(TAU(IA,2,ISP)+6.5_q,1._q)-0.5_q TAU(IA,3,ISP)=MOD(TAU(IA,3,ISP)+6.5_q,1._q)-0.5_q ! It could happen that (due to roundoff errors!) the surface coordinates ! are not uniquely -0.5 but probably +0.49999999999... . Ensure that all ! these coordinates are transformed to -0.5! IF ((ABS(TAU(IA,1,ISP)+0.5_q)) C_1 9 --> C_3 17 --> D_4 25 --> C_6v * ! 2 --> S_2 10 --> S_6 18 --> C_4v 26 --> D_3h * ! 3 --> C_2 11 --> D_3 19 --> D_2d 27 --> D_6h * ! 4 --> C_1h 12 --> C_3v 20 --> D_4h 28 --> T * ! 5 --> C_2h 13 --> D_3d 21 --> C_6 29 --> T_h * ! 6 --> D_2 14 --> C_4 22 --> C_3h 30 --> O * ! 7 --> C_2v 15 --> S_4 23 --> C_6h 31 --> T_d * ! 8 --> D_2h 16 --> C_4h 24 --> D_6 32 --> O_h * ! SFNAME is the explicit name in form of a string (Schoenflies). * ! * ! * !*********************************************************************** INTEGER S(3,3,48),PGIND,TRACE,DET CHARACTER*(*) SFNAME CHARACTER*4 GNAME(32) SAVE GNAME DATA GNAME /'C_1 ','S_2 ','C_2 ','C_1h','C_2h','D_2 ','C_2v', & & 'D_2h','C_3 ','S_6 ','D_3 ','C_3v','D_3d','C_4 ', & & 'S_4 ','C_4h','D_4 ','C_4v','D_2d','D_4h','C_6 ', & & 'C_3h','C_6h','D_6 ','C_6v','D_3h','D_6h','T ', & & 'T_h ','O ','T_d ','O_h '/ ! Need at least four characters to store the group name ... : IF (LEN(SFNAME)<4) CALL ERROR(' PGROUP', & & ' Variable SFNAME declared too short in the calling program!', & & LEN(SFNAME)) ! Trivial case: Only one symmetry operation (can only be E --> PGIND=1): IF (NROT==1) THEN PGIND=1 SFNAME=GNAME(PGIND) RETURN END IF ! There is an other trivial case: group O_h (PGIND=32), because it is ! the only group having 48 elements ... ! IF (NROT==48) THEN PGIND=32 SFNAME=GNAME(PGIND) RETURN END IF ! An other trivial case is group D_4h (PGIND=20), because it is the only ! group having 16 elements ... ! IF (NROT==16) THEN PGIND=20 SFNAME=GNAME(PGIND) RETURN END IF ! And finally there is a fourth trivial case: it is group C_3 (PGIND=9), ! because it is the only group having 3 elements ... ! IF (NROT==3) THEN PGIND=9 SFNAME=GNAME(PGIND) RETURN END IF ! All other groups need further investigations and detailed analysis ... ! First determine the type of elements and count them. Possible elements ! are E, I, C_2,3,4,6 and S_2,3,4,6 (S_2 = m), E is trivial and always ! present. The type of a symmetry operation can be identified simply by ! calculating the trace and the determinant of the rotation matrix. The ! combination of these two quantities is specific for specific elements: ! ! Element: E I C_2 C_3 C_4 C_6 S_2 S_6 S_4 S_3 ! Trace: +3 -3 -1 0 +1 +2 +1 0 -1 -2 ! Determinant: +1 -1 +1 +1 +1 +1 -1 -1 -1 -1 INVERS=0 NC2=0 NC3=0 NC4=0 NC6=0 NS2=0 NS6=0 NS4=0 NS3=0 DO 1 IR=1,NROT TRACE=S(1,1,IR)+S(2,2,IR)+S(3,3,IR) ! Found unity operator (trivial): IF (TRACE==3) GOTO 1 ! Found inversion ... IF (TRACE==(-3)) THEN INVERS=1 GOTO 1 END IF DET=S(1,1,IR)*(S(2,2,IR)*S(3,3,IR)-S(2,3,IR)*S(3,2,IR))+ & & S(1,2,IR)*(S(2,3,IR)*S(3,1,IR)-S(2,1,IR)*S(3,3,IR))+ & & S(1,3,IR)*(S(2,1,IR)*S(3,2,IR)-S(2,2,IR)*S(3,1,IR)) ! Found C_2: IF ((TRACE==(-1)).AND.(DET==1)) NC2=NC2+1 ! Found S_2: IF ((TRACE==1).AND.(DET==(-1))) NS2=NS2+1 ! Found C_3: IF ((TRACE==0).AND.(DET==1)) NC3=NC3+1 ! Found S_6: IF ((TRACE==0).AND.(DET==(-1))) NS6=NS6+1 ! Found C_4: IF ((TRACE==1).AND.(DET==1)) NC4=NC4+1 ! Found S_4: IF ((TRACE==(-1)).AND.(DET==(-1))) NS4=NS4+1 ! Found C_6: IF ((TRACE==2).AND.(DET==1)) NC6=NC6+1 ! Found S_3: IF ((TRACE==(-2)).AND.(DET==(-1))) NS3=NS3+1 1 CONTINUE ! Now we know which elements we have and so we know the group ... : IF (NROT==2) THEN ! Groups with 2 elements: IF (INVERS==1) THEN ! Contains inversion --> S_2 (PGIND=2): PGIND=2 SFNAME=GNAME(PGIND) RETURN END IF IF (NC2==1) THEN ! Contains twofold rotation --> C_2 (PGIND=3): PGIND=3 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==1) THEN ! Contains mirror plane --> C_1h (PGIND=4): PGIND=4 SFNAME=GNAME(PGIND) RETURN END IF END IF IF (NROT==4) THEN ! Groups with 4 elements: IF (INVERS==1) THEN ! Contains inversion --> C_2h (PGIND=5): PGIND=5 SFNAME=GNAME(PGIND) RETURN END IF IF (NC2==3) THEN ! Contains three twofold rotations --> D_2 (PGIND=6): PGIND=6 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==2) THEN ! Contains two mirror planes --> C_2v (PGIND=7): PGIND=7 SFNAME=GNAME(PGIND) RETURN END IF IF (NC4==1) THEN ! Contains fourfold rotation --> C_4 (PGIND=14): PGIND=14 SFNAME=GNAME(PGIND) RETURN END IF IF (NS4==1) THEN ! Contains fourfold improper rotation --> S_4 (PGIND=15): PGIND=15 SFNAME=GNAME(PGIND) RETURN END IF END IF IF (NROT==6) THEN ! Groups with 6 elements: IF (INVERS==1) THEN ! Contains inversion --> S_6 (PGIND=10): PGIND=10 SFNAME=GNAME(PGIND) RETURN END IF IF (NC2==3) THEN ! Contains three twofold rotations --> D_3 (PGIND=11): PGIND=11 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==3) THEN ! Contains three mirror planes --> C_3v (PGIND=12): PGIND=12 SFNAME=GNAME(PGIND) RETURN END IF IF (NC2==1) THEN ! Contains only one twofold rotations --> C_6 (PGIND=21): PGIND=21 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==1) THEN ! Contains only one mirror plane --> C_3h (PGIND=22): PGIND=22 SFNAME=GNAME(PGIND) RETURN END IF END IF IF (NROT==8) THEN ! Groups with 8 elements: IF (NS2==3) THEN ! Contains three mirror planes --> D_2h (PGIND=8): PGIND=8 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==1) THEN ! Contains one mirror planes --> C_4h (PGIND=16): PGIND=16 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==0) THEN ! Contains no mirror planes --> D_4 (PGIND=17): PGIND=17 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==4) THEN ! Contains four mirror planes --> C_4v (PGIND=18): PGIND=18 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==2) THEN ! Contains two mirror planes --> D_2d (PGIND=19): PGIND=19 SFNAME=GNAME(PGIND) RETURN END IF END IF IF (NROT==12) THEN ! Groups with 12 elements: IF (NS2==3) THEN ! Contains three mirror planes --> D_3d (PGIND=13): PGIND=13 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==1) THEN ! Contains one mirror planes --> C_6h (PGIND=23): PGIND=23 SFNAME=GNAME(PGIND) RETURN END IF IF (NC2==7) THEN ! Contains seven twofold rotations --> D_6 (PGIND=24): PGIND=24 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==6) THEN ! Contains six mirror planes --> C_6v (PGIND=25): PGIND=25 SFNAME=GNAME(PGIND) RETURN END IF IF (NS2==4) THEN ! Contains four mirror planes --> D_3h (PGIND=26): PGIND=26 SFNAME=GNAME(PGIND) RETURN END IF IF (NC3==8) THEN ! Contains eight threefold rotations --> T (PGIND=28): PGIND=28 SFNAME=GNAME(PGIND) RETURN END IF END IF IF (NROT==24) THEN ! Groups with 24 elements: IF (NC6==2) THEN ! Contains two sixfold rotations --> D_6h (PGIND=27): PGIND=27 SFNAME=GNAME(PGIND) RETURN END IF IF (INVERS==1) THEN ! Contains inversion --> T_h (PGIND=29): PGIND=29 SFNAME=GNAME(PGIND) RETURN END IF IF (NC4==6) THEN ! Contains six fourfold rotations --> O (PGIND=30): PGIND=30 SFNAME=GNAME(PGIND) RETURN END IF IF (NS4==6) THEN ! Contains six fourfold improper rotations --> T_d (PGIND=31): PGIND=31 SFNAME=GNAME(PGIND) RETURN END IF END IF ! Ready! RETURN END !*********************************************************************** ! * SUBROUTINE SGRGEN(SGEN,NGEN,S,NTOT) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine SGRGEN generates all point group symmetry operations from * ! the generation group. This routine is set up for INTEGER matrices. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! SGEN(9,NGEN) contains the generator matrices. * ! NGEN contains the number of generators. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! SGEN(9,48) contains the generated symmetry elements. * ! NTOT contains the number of generated symmetry operations. * ! * ! * !*********************************************************************** INTEGER SGEN(9,NGEN),S(9,48),H(9),HH(9),E(9),SIG(9) LOGICAL SGREQL SAVE E DATA E /1,0,0,0,1,0,0,0,1/ ! There is always one trivial element: the idendity ... CALL SGRCOP(E,S(1,1)) NTOT=1 ! Take all generators ... : DO 80 IGEN=1,NGEN CALL SGRCOP(SGEN(1,IGEN),SIG) ! Extend the group by all products with SIG(i)=SGEN(i,IGEN): DO 9 IG=1,NTOT IF (SGREQL(S(1,IG),SIG)) GOTO 80 9 CONTINUE ! Determine the order of the operation: CALL SGRCOP(SIG,H) DO 1 ITRY=1,100 IORD=ITRY IF (SGREQL(H,E)) GOTO 2 CALL SGRPRD(SIG,H,H) 1 CONTINUE ! Products of type 'G1*(SIG**P)*G2': 2 NNOW=NTOT DO 8 J=1,NTOT CALL SGRCOP(S(1,J),H) DO 10 IP=1,IORD-1 ! H=SIG**IP CALL SGRPRD(SIG,H,H) DO 11 I=1,NTOT ! HH=SYMOPS_I*(SIG**IP) CALL SGRPRD(S(1,I),H,HH) DO 12 K=1,NNOW IF (SGREQL(S(1,K),HH)) GOTO 11 12 CONTINUE 13 NNOW=NNOW+1 IF (NNOW>48) GOTO 99 CALL SGRCOP(HH,S(1,NNOW)) 11 CONTINUE 10 CONTINUE IF (J==1) N2=NNOW 8 CONTINUE ! Products with more than one sandwiched SIGMA-factor: M1=NTOT+1 M2=NNOW DO 20 I=2,50 DO 21 N=NTOT+1,N2 DO 21 M=M1,M2 CALL SGRPRD(S(1,N),S(1,M),H) DO 22 K=1,NNOW IF (SGREQL(S(1,K),H)) GOTO 21 22 CONTINUE NNOW=NNOW+1 IF (NNOW>48) GOTO 99 CALL SGRCOP(H,S(1,NNOW)) 21 CONTINUE IF (M2==NNOW) GOTO 25 M1=M2+1 M2=NNOW 20 CONTINUE 25 CONTINUE NTOT=NNOW 80 CONTINUE RETURN ! Error ! 99 CALL ERROR(' SGRGEN',' Too many elements',NNOW) RETURN END !*********************************************************************** ! * SUBROUTINE SGRPRD(S1,S2,S) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine SGRPRD multiplies two (INTEGER-) symmetry operators. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! S1(3,3) and S2(3,3) are the two matrices to be multiplied. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! S(3,3) contains the product matrix S1*S2. * ! * ! * !*********************************************************************** INTEGER S1(3,3),S2(3,3),S(3,3),SUM,HELP(3,3) DO 10 I=1,3 DO 10 J=1,3 SUM=0 DO 20 K=1,3 20 SUM=SUM+S1(I,K)*S2(K,J) 10 HELP(I,J)=SUM CALL SGRCOP(HELP,S) RETURN END !*********************************************************************** ! * SUBROUTINE SGRCOP(S1,S2) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine SGRPRD copies an (INTEGER-) symmetry operator. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! S1(9) is the matrix to be copied. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! S2(9) is the result matrix (copy of S1). * ! * ! * !*********************************************************************** INTEGER S1(9),S2(9) DO 10 I=1,9 S2(I)=S1(I) 10 CONTINUE RETURN END !*********************************************************************** ! * LOGICAL FUNCTION SGREQL(S1,S2) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Function SGREQL checks whether two (INTEGER-) symmetry operators * ! are identical or not ... . * ! * ! * ! Input parameters: * ! ----------------- * ! * ! S1(9) and S2(9) are two matrices to be compared. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! SGREQL is the result (TRUE if S1=S2 and FALSE else). * ! * ! * !*********************************************************************** INTEGER S1(9),S2(9) SGREQL=.FALSE. DO 10 I=1,9 IF (S1(I)/=S2(I)) RETURN 10 CONTINUE SGREQL=.TRUE. RETURN END !*********************************************************************** ! * SUBROUTINE INVGRP(S,NROT,INVMAP) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine INVGRP searches the inverse operators of a given set of * ! 'rotation' matrices and creates a map where to find the inverse * ! operator ('storage index' of the inverse operator). * ! * ! * ! Input parameters: * ! ----------------- * ! * ! S(3,3,48) contains the input matrices. * ! NROT gives the number of matrices in S. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! INVMAP(48) contains the index table for the inverse operators. * ! * ! * !*********************************************************************** INTEGER S(3,3,48),G(9),INVMAP(48),E(9) LOGICAL SGREQL SAVE E DATA E /1,0,0,0,1,0,0,0,1/ ! For all symmetry operations ... : DO 200 IROT=1,NROT ! Test all symmetry operations to find the inverse ... : DO 100 JROT=1,NROT ! Multiply the 'rotation' matrix with the test 'rotation' matrix ... CALL SGRPRD(S(1,1,IROT),S(1,1,JROT),G(1)) ! ... and if the result is equal to the unit operator (unit matrix) ! then we found the inverse operator ... . IF (SGREQL(G(1),E(1))) THEN ! Store the result ... INVMAP(IROT)=JROT ! ... and leave the 'test loop' ... . GOTO 150 END IF 100 CONTINUE 150 CONTINUE 200 CONTINUE RETURN END !*********************************************************************** ! * SUBROUTINE SGRCON(SA,SB,NROT,A1,A2,A3,B1,B2,B3) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine SGRCON converts (INTEGER-) rotation matrices given in some * ! arbitrary basis A to its representation in some other basis B. * ! WARNING: There is one restriction: the two basis sets A and B must * ! be 'compatible' (same crystal class, inverse basis sets, * ! or B = integer linear combination of A) in order to get * ! again INTEGER output matrices!! Otherwise use SGRTRF! * ! * ! * ! Input parameters: * ! ----------------- * ! * ! SA(3,3,48) contains the input matrices (basis A1,A2,A3). * ! NROT gives the number of matrices in SA (and hence in SB). * ! * ! A1(3),A2(3),A3(3) contain the first set of basis vectors. * ! B1(3),B2(3),B3(3) contain the second set of basis vectors. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! SB(3,3,48) contains the converted matrices (basis B1,B2,B3). * ! * ! * !*********************************************************************** INTEGER SA(3,3,48),SB(3,3,48) DIMENSION ADOTBI(3,3),BDOTAI(3,3),PROD(3,3) DIMENSION A1(3),A2(3),A3(3),B1(3),B2(3),B3(3),R(3,3) DIMENSION AI1(3),AI2(3),AI3(3),BI1(3),BI2(3),BI3(3) SAVE TINY DATA TINY / MY_TINY / ! Calculate reciprocal basis vectors for both sets of direct vectors: CALL RECIPS(1._q,A1(1),A2(1),A3(1),AI1(1),AI2(1),AI3(1)) CALL RECIPS(1._q,B1(1),B2(1),B3(1),BI1(1),BI2(1),BI3(1)) ! Set up the matrices of scalar products (a_i*bi_j) and (b_i*ai_j): BDOTAI(1,1)=AI1(1)*B1(1)+AI1(2)*B1(2)+AI1(3)*B1(3) BDOTAI(2,2)=AI2(1)*B2(1)+AI2(2)*B2(2)+AI2(3)*B2(3) BDOTAI(3,3)=AI3(1)*B3(1)+AI3(2)*B3(2)+AI3(3)*B3(3) BDOTAI(2,1)=AI1(1)*B2(1)+AI1(2)*B2(2)+AI1(3)*B2(3) BDOTAI(3,1)=AI1(1)*B3(1)+AI1(2)*B3(2)+AI1(3)*B3(3) BDOTAI(1,2)=AI2(1)*B1(1)+AI2(2)*B1(2)+AI2(3)*B1(3) BDOTAI(3,2)=AI2(1)*B3(1)+AI2(2)*B3(2)+AI2(3)*B3(3) BDOTAI(1,3)=AI3(1)*B1(1)+AI3(2)*B1(2)+AI3(3)*B1(3) BDOTAI(2,3)=AI3(1)*B2(1)+AI3(2)*B2(2)+AI3(3)*B2(3) ADOTBI(1,1)=BI1(1)*A1(1)+BI1(2)*A1(2)+BI1(3)*A1(3) ADOTBI(2,2)=BI2(1)*A2(1)+BI2(2)*A2(2)+BI2(3)*A2(3) ADOTBI(3,3)=BI3(1)*A3(1)+BI3(2)*A3(2)+BI3(3)*A3(3) ADOTBI(2,1)=BI1(1)*A2(1)+BI1(2)*A2(2)+BI1(3)*A2(3) ADOTBI(3,1)=BI1(1)*A3(1)+BI1(2)*A3(2)+BI1(3)*A3(3) ADOTBI(1,2)=BI2(1)*A1(1)+BI2(2)*A1(2)+BI2(3)*A1(3) ADOTBI(3,2)=BI2(1)*A3(1)+BI2(2)*A3(2)+BI2(3)*A3(3) ADOTBI(1,3)=BI3(1)*A1(1)+BI3(2)*A1(2)+BI3(3)*A1(3) ADOTBI(2,3)=BI3(1)*A2(1)+BI3(2)*A2(2)+BI3(3)*A2(3) ! For all (rotation-) matrices: DO 100 I=1,NROT R(1,1)=FLOAT(SA(1,1,I)) R(2,1)=FLOAT(SA(2,1,I)) R(3,1)=FLOAT(SA(3,1,I)) R(1,2)=FLOAT(SA(1,2,I)) R(2,2)=FLOAT(SA(2,2,I)) R(3,2)=FLOAT(SA(3,2,I)) R(1,3)=FLOAT(SA(1,3,I)) R(2,3)=FLOAT(SA(2,3,I)) R(3,3)=FLOAT(SA(3,3,I)) ! Perform the transformation (matrix multiplications): DO 101 II=1,3 DO 101 JJ=1,3 SUM=0._q DO 201 KK=1,3 201 SUM=SUM+BDOTAI(II,KK)*R(KK,JJ) 101 PROD(II,JJ)=SUM DO 102 II=1,3 DO 102 JJ=1,3 SUM=0._q DO 202 KK=1,3 202 SUM=SUM+PROD(II,KK)*ADOTBI(KK,JJ) 102 R(II,JJ)=SUM ! Ready! Copy the result into SB ... : SB(1,1,I)=NINT(R(1,1)) SB(2,1,I)=NINT(R(2,1)) SB(3,1,I)=NINT(R(3,1)) SB(1,2,I)=NINT(R(1,2)) SB(2,2,I)=NINT(R(2,2)) SB(3,2,I)=NINT(R(3,2)) SB(1,3,I)=NINT(R(1,3)) SB(2,3,I)=NINT(R(2,3)) SB(3,3,I)=NINT(R(3,3)) ! Error check: We must end up with INTEGER matrices ... ! IF (ABS(R(1,1)-SB(1,1,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(2,1)-SB(2,1,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(3,1)-SB(3,1,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(1,2)-SB(1,2,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(2,2)-SB(2,2,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(3,2)-SB(3,2,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(1,3)-SB(1,3,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(2,3)-SB(2,3,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) IF (ABS(R(3,3)-SB(3,3,I))>TINY) CALL ERROR(' SGRCON', & & ' Found some non-integer element in rotation matrix',I) 100 CONTINUE ! Good bye ... : RETURN END !*********************************************************************** ! * SUBROUTINE SGRTRF(SA,SB,NROT,A1,A2,A3,B1,B2,B3) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine SGRTRF converts (REAL-) rotation matrices given in some * ! arbitrary basis A to its representation in some other basis B. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! SA(3,3,48) contains the input matrices (basis A1,A2,A3). * ! NROT gives the number of matrices in SA (and hence in SB). * ! * ! A1(3),A2(3),A3(3) contain the first set of basis vectors. * ! B1(3),B2(3),B3(3) contain the second set of basis vectors. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! SB(3,3,48) contains the converted matrices (basis B1,B2,B3). * ! * ! * !*********************************************************************** DIMENSION SA(3,3,48),SB(3,3,48) DIMENSION ADOTBI(3,3),BDOTAI(3,3),PROD(3,3) DIMENSION A1(3),A2(3),A3(3),B1(3),B2(3),B3(3),R(3,3) DIMENSION AI1(3),AI2(3),AI3(3),BI1(3),BI2(3),BI3(3) ! Calculate reciprocal basis vectors for both sets of direct vectors: CALL RECIPS(1._q,A1(1),A2(1),A3(1),AI1(1),AI2(1),AI3(1)) CALL RECIPS(1._q,B1(1),B2(1),B3(1),BI1(1),BI2(1),BI3(1)) ! Set up the matrices of scalar products (a_i*bi_j) and (b_i*ai_j): BDOTAI(1,1)=AI1(1)*B1(1)+AI1(2)*B1(2)+AI1(3)*B1(3) BDOTAI(2,2)=AI2(1)*B2(1)+AI2(2)*B2(2)+AI2(3)*B2(3) BDOTAI(3,3)=AI3(1)*B3(1)+AI3(2)*B3(2)+AI3(3)*B3(3) BDOTAI(2,1)=AI1(1)*B2(1)+AI1(2)*B2(2)+AI1(3)*B2(3) BDOTAI(3,1)=AI1(1)*B3(1)+AI1(2)*B3(2)+AI1(3)*B3(3) BDOTAI(1,2)=AI2(1)*B1(1)+AI2(2)*B1(2)+AI2(3)*B1(3) BDOTAI(3,2)=AI2(1)*B3(1)+AI2(2)*B3(2)+AI2(3)*B3(3) BDOTAI(1,3)=AI3(1)*B1(1)+AI3(2)*B1(2)+AI3(3)*B1(3) BDOTAI(2,3)=AI3(1)*B2(1)+AI3(2)*B2(2)+AI3(3)*B2(3) ADOTBI(1,1)=BI1(1)*A1(1)+BI1(2)*A1(2)+BI1(3)*A1(3) ADOTBI(2,2)=BI2(1)*A2(1)+BI2(2)*A2(2)+BI2(3)*A2(3) ADOTBI(3,3)=BI3(1)*A3(1)+BI3(2)*A3(2)+BI3(3)*A3(3) ADOTBI(2,1)=BI1(1)*A2(1)+BI1(2)*A2(2)+BI1(3)*A2(3) ADOTBI(3,1)=BI1(1)*A3(1)+BI1(2)*A3(2)+BI1(3)*A3(3) ADOTBI(1,2)=BI2(1)*A1(1)+BI2(2)*A1(2)+BI2(3)*A1(3) ADOTBI(3,2)=BI2(1)*A3(1)+BI2(2)*A3(2)+BI2(3)*A3(3) ADOTBI(1,3)=BI3(1)*A1(1)+BI3(2)*A1(2)+BI3(3)*A1(3) ADOTBI(2,3)=BI3(1)*A2(1)+BI3(2)*A2(2)+BI3(3)*A2(3) ! For all (rotation-) matrices: DO 100 I=1,NROT R(1,1)=SA(1,1,I) R(2,1)=SA(2,1,I) R(3,1)=SA(3,1,I) R(1,2)=SA(1,2,I) R(2,2)=SA(2,2,I) R(3,2)=SA(3,2,I) R(1,3)=SA(1,3,I) R(2,3)=SA(2,3,I) R(3,3)=SA(3,3,I) ! Perform the transformation (matrix multiplications): DO 101 II=1,3 DO 101 JJ=1,3 SUM=0._q DO 201 KK=1,3 201 SUM=SUM+BDOTAI(II,KK)*R(KK,JJ) 101 PROD(II,JJ)=SUM DO 102 II=1,3 DO 102 JJ=1,3 SUM=0._q DO 202 KK=1,3 202 SUM=SUM+PROD(II,KK)*ADOTBI(KK,JJ) 102 R(II,JJ)=SUM ! Ready! Copy the result into SB ... : SB(1,1,I)=R(1,1) SB(2,1,I)=R(2,1) SB(3,1,I)=R(3,1) SB(1,2,I)=R(1,2) SB(2,2,I)=R(2,2) SB(3,2,I)=R(3,2) SB(1,3,I)=R(1,3) SB(2,3,I)=R(2,3) SB(3,3,I)=R(3,3) 100 CONTINUE ! Good bye ... : RETURN END !*********************************************************************** ! * SUBROUTINE VECCON(VA,VB,NVEC,A1,A2,A3,B1,B2,B3) USE prec IMPLICIT REAL(q) (A-H,O-Z) ! * !*********************************************************************** ! * ! Routine VECCON converts a set of vectors represented in some * ! arbitrary basis A to its representation in some other basis B. * ! * ! * ! Input parameters: * ! ----------------- * ! * ! VA(3,NVEC) contains the input vectors (basis A1,A2,A3). * ! NVEC gives the number of vectors in VA (and hence in VB). * ! * ! A1(3),A2(3),A3(3) contain the first set of basis vectors. * ! B1(3),B2(3),B3(3) contain the second set of basis vectors. * ! * ! * ! Output parameters: * ! ------------------ * ! * ! VB(3,NVEC) contains the converted vectors (basis B1,B2,B3). * ! * ! * !*********************************************************************** DIMENSION VA(3,NVEC),VB(3,NVEC) DIMENSION A1(3),A2(3),A3(3),B1(3),B2(3),B3(3) DIMENSION BI1(3),BI2(3),BI3(3),V(3),R(3) ! Calculate reciprocal basis vectors for second set of direct vectors: CALL RECIPS(1._q,B1(1),B2(1),B3(1),BI1(1),BI2(1),BI3(1)) ! For all vectors ... : DO 100 I=1,NVEC ! Copy the vector: V(1)=VA(1,I) V(2)=VA(2,I) V(3)=VA(3,I) ! Perform the transformation: R(1)=A1(1)*V(1)+A2(1)*V(2)+A3(1)*V(3) R(2)=A1(2)*V(1)+A2(2)*V(2)+A3(2)*V(3) R(3)=A1(3)*V(1)+A2(3)*V(2)+A3(3)*V(3) V(1)=BI1(1)*R(1)+BI1(2)*R(2)+BI1(3)*R(3) V(2)=BI2(1)*R(1)+BI2(2)*R(2)+BI2(3)*R(3) V(3)=BI3(1)*R(1)+BI3(2)*R(2)+BI3(3)*R(3) ! Store the result: VB(1,I)=V(1) VB(2,I)=V(2) VB(3,I)=V(3) 100 CONTINUE ! Good bye ... : RETURN END