TABLE OF CONTENTS
- ABINIT/m_symfind
- m_symfind/symanal
- m_symfind/symbrav
- m_symfind/symfind
- m_symfind/symfind_expert
- m_symfind/symlatt
- m_symfind/symspgr
ABINIT/m_symfind [ Modules ]
NAME
m_symfind
FUNCTION
Symmetry finder high-level API.
COPYRIGHT
Copyright (C) 2000-2024 ABINIT group (XG, RC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_symfind 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_symlist 28 29 use m_symtk, only : chkprimit, mati3inv, matr3inv, symrelrot, symdet, symcharac, holocell, & 30 smallprim, print_symmetries, sg_multable, symatm, symmetrize_tnons, symmetrize_xred 31 use m_geometry, only : acrossb, xred2xcart 32 use m_spgdata, only : getptgroupma, symptgroup, spgdata 33 34 implicit none 35 36 private
m_symfind/symanal [ Functions ]
[ Top ] [ m_symfind ] [ Functions ]
NAME
symanal
FUNCTION
Find the space group, Bravais lattice, including Shubnikov characteristics from the list of symmetries (including magnetic characteristics), and lattice parameters Warning: the recognition of the space group might not yet work for the Shubnikov group of type IV
INPUTS
chkprim= if 1 then stop if the cell is not primitive msym=default maximal number of symmetries nsym=actual number of symmetries rprimd(3,3)=dimensional primitive translations for real space (bohr) symafm(1:msym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,1:msym)=symmetry operations in real space in terms of primitive translations tnons(3,1:msym)=nonsymmorphic translations for symmetry operations tolsym=tolerance for the symmetry operations verbose= if true, will list the symmetry operation labels
OUTPUT
bravais(11)=characteristics of Bravais lattice (see symlatt.F90) genafm(3)=magnetic translation generator (in case of Shubnikov group type IV) ptgroupma = magnetic point group number spgroup=symmetry space group
SOURCE
885 subroutine symanal(bravais,chkprim,genafm,msym,nsym,ptgroupma,rprimd,spgroup,symafm,symrel,tnons,tolsym,verbose) 886 887 !Arguments ------------------------------------ 888 !scalars 889 integer,intent(in) :: chkprim,msym,nsym 890 integer,intent(out) :: ptgroupma,spgroup 891 real(dp),intent(in) :: tolsym 892 logical,optional,intent(in) :: verbose 893 !arrays 894 integer,intent(out) :: bravais(11) 895 integer,intent(in) :: symafm(msym),symrel(3,3,msym) 896 real(dp),intent(in) :: rprimd(3,3) 897 real(dp),intent(in) :: tnons(3,msym) 898 real(dp),intent(out) :: genafm(3) 899 900 !Local variables------------------------------- 901 !scalars 902 integer, parameter :: maxsym=192 903 ! In this routine, maxsym is used either to determine the ptsymrel from the rprimd (routine symlatt), 904 ! so, it might be up to 192 = 4*48 for FCC, and also to define the maximum number of symmetry operation labels, 905 ! but only in case the cell is primitive, which gives the same upper bound. Thus in this routine, msym might 906 ! be equal to nsym. 907 integer :: iholohedry_nomagn,isym,isym_nomagn,multi 908 integer :: nptsym,nsym_nomagn,shubnikov 909 logical :: verbose_ 910 character(len=5) :: ptgroup,ptgroupha 911 character(len=500) :: msg 912 !arrays 913 integer :: identity(3,3) 914 integer,allocatable :: ptsymrel(:,:,:),symrel_nomagn(:,:,:) 915 real(dp),allocatable :: tnons_nomagn(:,:) 916 character(len=128) :: labels(maxsym) 917 918 ! ************************************************************************* 919 920 !DEBUG 921 !write(std_out,*)' symanal : enter' 922 !write(std_out,*)' symanal : chkprim =',chkprim 923 !write(std_out,*)' symanal : nsym=',nsym 924 !do isym=1,nsym 925 ! write(std_out,*)' symanal : symrel=',symrel(1:3,1:3,isym) 926 !enddo 927 !ENDDEBUG 928 929 verbose_=.false. 930 if(present(verbose))then 931 verbose_=verbose 932 endif 933 934 !This routine finds the Bravais characteristics, without actually 935 !looking at the symmetry operations. 936 ABI_MALLOC(ptsymrel,(3,3,maxsym)) 937 call symlatt(bravais,dev_null,maxsym,nptsym,ptsymrel,rprimd,tolsym) 938 ABI_FREE(ptsymrel) 939 940 !Check whether the cell is primitive or not. 941 call chkprimit(chkprim,multi,nsym,symafm,symrel) 942 943 spgroup=0 ; ptgroupma=0 ; genafm(:)=zero 944 945 if(multi>1)then ! Modify bravais if the cell is not primitive ; no determination of the space group 946 bravais(1)=-bravais(1) 947 else 948 949 ! The cell is primitive, so that the space group can be 950 ! determined. Need to distinguish Fedorov and Shubnikov groups. 951 ! Do not distinguish Shubnikov types I and II. 952 ! Also identify genafm, in case of Shubnikov type IV 953 identity(:,:)=reshape((/1,0,0,0,1,0,0,0,1/),(/3,3/)) 954 shubnikov=1 955 do isym=1,nsym 956 if(symafm(isym)==-1)then 957 shubnikov=3 958 if(sum(abs(symrel(:,:,isym)-identity(:,:)))==0)then 959 shubnikov=4 960 genafm(:)=tnons(:,isym) 961 ! DEBUG 962 ! write(std_out,*)' isym=',isym 963 ! write(std_out,*)' symrel(:,:,isym)',symrel(:,:,isym) 964 ! write(std_out,*)' tnons(:,isym)',tnons(:,isym) 965 ! write(std_out,*)' symafm(isym)',symafm(isym) 966 ! ENDDEBUG 967 exit 968 end if 969 end if 970 end do 971 972 if(shubnikov/=1)then 973 if(shubnikov==3)write(msg, '(a)' )' Shubnikov space group type III' 974 if(shubnikov==4)write(msg, '(a)' )' Shubnikov space group type IV' 975 call wrtout(std_out,msg) 976 end if 977 978 if(shubnikov==1 .or. shubnikov==3)then 979 ! Find the correct Bravais characteristics and point group 980 ! Should also be used for Shubnikov groups of type IV ... 981 call symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym) 982 983 ! Find the space group 984 call symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym) 985 986 if(verbose_)then 987 do isym=1,nsym 988 write(msg,'(a,i3,2a)')' symanal : the symmetry operation no. ',isym,' is ',trim(labels(isym)) 989 call wrtout(std_out,msg) 990 enddo 991 endif 992 993 end if 994 995 if(shubnikov/=1)then 996 997 ! Determine nonmagnetic symmetry operations 998 nsym_nomagn=nsym/2 999 ABI_MALLOC(symrel_nomagn,(3,3,nsym_nomagn)) 1000 ABI_MALLOC(tnons_nomagn,(3,nsym_nomagn)) 1001 isym_nomagn=0 1002 do isym=1,nsym 1003 if(symafm(isym)==1)then 1004 isym_nomagn=isym_nomagn+1 1005 symrel_nomagn(:,:,isym_nomagn)=symrel(:,:,isym) 1006 tnons_nomagn(:,isym_nomagn)=tnons(:,isym) 1007 end if 1008 end do 1009 1010 if(shubnikov==3)then 1011 1012 ! DEBUG 1013 ! write(std_out,*)' symanal : will enter symbrav with halved symmetry set' 1014 ! write(std_out,*)' Describe the different symmetry operations (index,symrel,tnons,symafm)' 1015 ! do isym=1,nsym_nomagn 1016 ! write(std_out,'(i3,2x,9i3,3es12.2,i3)')isym,symrel_nomagn(:,:,isym),tnons_nomagn(:,isym) 1017 ! end do 1018 ! ENDDEBUG 1019 1020 ! Find the point group of the halved symmetry set 1021 call symptgroup(iholohedry_nomagn,nsym_nomagn,ptgroupha,symrel_nomagn) 1022 1023 ! Deduce the magnetic point group (ptgroupma) from ptgroup and ptgroupha 1024 call getptgroupma(ptgroup,ptgroupha,ptgroupma) 1025 1026 else if(shubnikov==4)then 1027 1028 ! Find the Fedorov space group of the halved symmetry set 1029 call symspgr(bravais,labels,nsym_nomagn,spgroup,symrel_nomagn,tnons_nomagn,tolsym) 1030 1031 ! The magnetic translation generator genafm has already been determined 1032 ! write(std_out,*)' genafm =',genafm, ' spgroup=',spgroup 1033 1034 if(verbose_)then 1035 1036 write(msg, '(a)' )' Select only the non-magnetic symmetry operations ' 1037 call wrtout(std_out,msg) 1038 1039 do isym=1,nsym 1040 if(symafm(isym)==1)then 1041 isym_nomagn=isym_nomagn+1 1042 write(msg,'(a,i3,2a)')' symspgr : the symmetry operation no. ',isym,' is ',trim(labels(isym_nomagn)) 1043 call wrtout(std_out,msg) 1044 endif 1045 enddo 1046 endif 1047 1048 end if 1049 1050 ABI_FREE(symrel_nomagn) 1051 ABI_FREE(tnons_nomagn) 1052 end if ! Shubnikov groups 1053 1054 end if 1055 1056 !DEBUG 1057 !write(std_out,'(a)') ' symanal : exit ' 1058 !ENDDEBUG 1059 1060 end subroutine symanal
m_symfind/symbrav [ Functions ]
[ Top ] [ m_symfind ] [ Functions ]
NAME
symbrav
FUNCTION
From the list of symmetry operations, and the lattice vectors, determine the Bravais information (including the holohedry, the centering, the coordinate of the primitive vectors in the conventional vectors), as well as the point group.
INPUTS
msym=dimension of symrel nsym=actual number of symmetries rprimd(3,3)=dimensional primitive translations for real space (bohr) symrel(3,3,msym)=symmetry operations in real space in terms of primitive translations tolsym=tolerance for the symmetries
OUTPUT
bravais(11): bravais(1)=iholohedry
bravais(2)=center
bravais(3:11)=coordinates of rprimd in the axes
of the conventional bravais lattice (*2 if center/=0)
ptgroup=symmetry point group
[axis(3)]=Invariant axis in the conventional vector coordinates
Set to (/0,0,0/) if the lattice belongs to the same holohedry as the lattice+atoms (+electric field + ...).
SOURCE
1092 subroutine symbrav(bravais,msym,nsym,ptgroup,rprimd,symrel,tolsym,axis) 1093 1094 !Arguments ------------------------------------ 1095 !scalars 1096 integer,intent(in) :: msym,nsym 1097 real(dp),intent(in) :: tolsym 1098 character(len=5),intent(out) :: ptgroup 1099 !arrays 1100 integer,intent(in) :: symrel(3,3,msym) 1101 integer,optional,intent(out) :: axis(3) 1102 integer,intent(out) :: bravais(11) 1103 real(dp),intent(in) :: rprimd(3,3) 1104 1105 !Local variables------------------------------- 1106 !scalars 1107 integer :: iaxis,ii,bravais1now,ideform,iholohedry,invariant,isym 1108 integer :: jaxis,next_stage,nptsym,problem,maxsym 1109 integer, parameter :: naxes_ortho=22, naxes_hexa=7 1110 real(dp) :: norm,scprod 1111 character(len=500) :: msg 1112 !arrays 1113 integer :: identity(3,3),axis_trial(3),hexa_axes(3,naxes_hexa),ortho_axes(3,naxes_ortho) 1114 integer,allocatable :: ptsymrel(:,:,:),symrelconv(:,:,:) 1115 real(dp) :: axes(3,3),axis_cart(3),axis_red(3) 1116 real(dp) :: rprimdconv(3,3),rprimdtry(3,3),rprimdnow(3,3) 1117 real(dp) :: rprimdconv_invt(3,3) 1118 1119 !************************************************************************** 1120 1121 !DEBUG 1122 !write(std_out,*)' symbrav : enter ' 1123 !call flush(std_out) 1124 !ENDDEBUG 1125 1126 identity(:,:)=0 1127 identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 1128 1129 ortho_axes(:,:)=0 1130 ortho_axes(1,1)=1 1131 ortho_axes(2,2)=1 1132 ortho_axes(3,3)=1 1133 ortho_axes(:,4)=(/0,1,1/) 1134 ortho_axes(:,5)=(/1,0,1/) 1135 ortho_axes(:,6)=(/1,1,0/) 1136 ortho_axes(:,7)=(/0,1,-1/) 1137 ortho_axes(:,8)=(/-1,0,1/) 1138 ortho_axes(:,9)=(/1,-1,0/) 1139 ortho_axes(:,10)=(/0,1,2/) 1140 ortho_axes(:,11)=(/2,0,1/) 1141 ortho_axes(:,12)=(/1,2,0/) 1142 ortho_axes(:,13)=(/1,1,1/) 1143 ortho_axes(:,14)=(/-1,1,1/) 1144 ortho_axes(:,15)=(/1,-1,1/) 1145 ortho_axes(:,16)=(/1,1,-1/) 1146 ortho_axes(:,17)=(/2,1,1/) 1147 ortho_axes(:,18)=(/1,2,1/) 1148 ortho_axes(:,19)=(/1,1,2/) 1149 ortho_axes(:,20)=(/2,1,-1/) 1150 ortho_axes(:,21)=(/-1,2,1/) 1151 ortho_axes(:,22)=(/1,-1,2/) 1152 1153 hexa_axes(:,:)=0 1154 hexa_axes(1,1)=1 1155 hexa_axes(2,2)=1 1156 hexa_axes(3,3)=1 1157 hexa_axes(:,4)=(/1,-1,0/) 1158 hexa_axes(:,5)=(/2,1,0/) 1159 hexa_axes(:,6)=(/1,1,0/) 1160 hexa_axes(:,7)=(/1,2,0/) 1161 1162 !Determine the point group from the list of symmetry operations. 1163 !Also determine the holohedry, up to one undeterminacy : hR versus hP 1164 call symptgroup(iholohedry,nsym,ptgroup,symrel) 1165 1166 !DEBUG 1167 !write(std_out,*)' symbrav, after symptgroup: nsym=',nsym 1168 !call flush(std_out) 1169 !write(std_out,*)' symbrav: symrel=' 1170 !do isym=1,nsym 1171 ! write(std_out,'(9i4)')symrel(:,:,isym) 1172 !enddo 1173 !write(std_out,*)' symbrav: iholohedry=',iholohedry 1174 !call flush(std_out) 1175 !ENDDEBUG 1176 1177 !Loop over trial deformations 1178 !This is needed in case the Bravais lattice determination from the lattice vectors 1179 !has a higher holohedry than the real one, in which the symmetry 1180 !operations for the atoms (or electric field, etc) are taken into account 1181 iaxis=0 1182 invariant=0 1183 next_stage=0 1184 rprimdnow(:,:)=rprimd(:,:) 1185 rprimdtry(:,:)=rprimd(:,:) 1186 ABI_MALLOC(symrelconv,(3,3,nsym)) 1187 1188 !At most will have to try naxes_ortho*5 deformations (naxes_ortho axes, five stages) 1189 !First, test whether the current recognition of Bravais lattice is problematic (iholohedry differs from bravais(1)). 1190 !Then, if there is a problem, test different deformations of rprimd, one after the other. 1191 !For each try, there is a new rprimdtry from the different set of deformation, that is generated later in the loop 1192 !Also bravais1now and rprimdnow might changed (and progressively lowered). 1193 !The latter change induces at most 5 stages for the computation (cubic->tetragonal->orthorhombic->monoclinic->triclinic). 1194 !After an upgrade of bravais1now and rprimdnow, one has to restart the full set of deformations. 1195 !Not sure that this procedure resolves all cases, but seems to work on >40000 inaccurate POSCAR files. 1196 do ideform=1,naxes_ortho*5 1197 1198 maxsym=max(192,msym) 1199 ABI_MALLOC(ptsymrel,(3,3,maxsym)) 1200 call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,tolsym) 1201 ABI_FREE(ptsymrel) 1202 1203 !DEBUG 1204 !write(std_out,*)' symbrav: inside loop with ideform,iaxis=',ideform,iaxis 1205 !write(std_out,'(a,9f12.6)')' rprimdtry=',rprimdtry(:,:) 1206 !write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2) 1207 !call flush(std_out) 1208 !ENDDEBUG 1209 1210 1211 ! Examine the agreement with bravais(1) 1212 ! Warning : might change Bravais lattice hR to hP, if hexagonal axes 1213 problem=0 1214 select case (bravais(1)) 1215 case(7) 1216 if(iholohedry<6)problem=1 1217 if(iholohedry==6)problem=2 1218 case(6) 1219 if(iholohedry<4)problem=1 1220 if(iholohedry==7 .or. iholohedry==4)problem=2 1221 ! Here, change hR into hP 1222 if(iholohedry==5)iholohedry=6 1223 case(5) 1224 if(iholohedry<4)problem=1 1225 if(iholohedry==7 .or. iholohedry==6 .or. iholohedry==4)problem=2 1226 case(4) 1227 if(iholohedry<4)problem=1 1228 if(iholohedry>4)problem=2 1229 case(3) 1230 if(iholohedry<3)problem=1 1231 if(iholohedry>3)problem=2 1232 case(2) 1233 if(iholohedry<2)problem=1 1234 if(iholohedry>2)problem=2 1235 case(1) 1236 if(iholohedry>1)problem=2 1237 end select 1238 1239 ! This is the usual situation, in which the lattice belong to the same holohedry 1240 ! as the lattice+atoms (+electric field + ...) 1241 if(problem==0)exit 1242 1243 if(problem==2)then 1244 if(iaxis==0)then 1245 write(msg, '(3a,i3,3a,i3,7a)' )& 1246 'The Bravais lattice determined only from the primitive',ch10,& 1247 'vectors (rprim or angdeg), bravais(1)=',bravais(1),', is not compatible',ch10,& 1248 'with the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,& 1249 'account the symmetry operations. This might be due to an insufficient',ch10,& 1250 'number of digits in the specification of rprim (at least 10),',ch10,& 1251 'or to an erroneous rprim or angdeg. If this is not the case, then ...' 1252 ABI_BUG(msg) 1253 end if 1254 if(iaxis==1)then 1255 write(msg, '(3a,3i3,2a,i3,2a,i3)' )& 1256 'Could not succeed to determine the bravais lattice',ch10,& 1257 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,& 1258 'bravais(1)=',bravais(1),ch10,& 1259 'iholohedry=',iholohedry 1260 ABI_BUG(msg) 1261 end if 1262 ! Try to increase tolsym to find the Bravais lattice. 1263 maxsym=max(192,msym) 1264 ABI_MALLOC(ptsymrel,(3,3,maxsym)) 1265 !DEBUG 1266 ! write(6,*)' symbrav : will call symlatt, 3*tolsym=',3*tolsym 1267 !ENDDEBUG 1268 call symlatt(bravais,std_out,maxsym,nptsym,ptsymrel,rprimdtry,3*tolsym) 1269 ABI_FREE(ptsymrel) 1270 if(bravais(1)==iholohedry)then 1271 ! Succeeded 1272 exit 1273 else 1274 write(msg, '(3a,3i3,2a,i3,2a,i3)' )& 1275 'Could not succeed to determine the bravais lattice, even after considering a larger tolsym',ch10,& 1276 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,& 1277 'bravais(1)=',bravais(1),ch10,& 1278 'iholohedry=',iholohedry 1279 ABI_BUG(msg) 1280 end if 1281 end if 1282 1283 if(problem==1)then ! One is left with the problem=1 case, basically iholohedry is lower than bravais(1) 1284 if(iaxis==0)then 1285 write(msg, '(a,a,a,i3,a,a,a,i3,a,a,a)' )& 1286 'The Bravais lattice determined only from the primitive',ch10,& 1287 'vectors, bravais(1)=',bravais(1),', is more symmetric',ch10,& 1288 'than the real one, iholohedry=',iholohedry,', obtained by taking into',ch10,& 1289 'account the atomic positions. Start deforming the primitive vector set.' 1290 ABI_COMMENT(msg) 1291 next_stage=1 1292 else if(iaxis/=0)then 1293 if(bravais(1)<bravais1now)then 1294 write(msg, '(3a,i3,3a,i3,2a)' )& 1295 'The Bravais lattice determined from modified primitive',ch10,& 1296 'vectors, bravais(1)=',bravais(1),', has a lower symmetry than before,',ch10,& 1297 'but is still more symmetric than the real one, iholohedry=',iholohedry,ch10,& 1298 'obtained by taking into account the atomic positions.' 1299 ABI_COMMENT(msg) 1300 next_stage=1 1301 else if(iaxis==1)then 1302 write(msg, '(3a,3i3,2a,i3,2a,i3)' )& 1303 'Could not succeed to determine the bravais lattice',ch10,& 1304 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,& 1305 'bravais(1)=',bravais(1),ch10,& 1306 'iholohedry=',iholohedry 1307 ABI_BUG(msg) 1308 end if 1309 end if 1310 end if ! problem==1 1311 1312 ! One is here when problem=1 (iholohedry < bravais(1)) and either 1313 ! - iaxis=0 (no deformation has been tried yet), 1314 ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) < bravais1now 1315 ! - some deformation iaxis has been tried giving bravais(1), but iholohedry < bravais(1) = bravais1now and iaxis/=1 . 1316 ! Also, note that next_stage is still 0 when bravais(1)=bravais1now . 1317 ! The loop has been ended (so, the search failed) when bravais(1)=bravais1now and iaxis==1. 1318 1319 if(next_stage==1)then 1320 bravais1now=bravais(1) 1321 rprimdnow(:,:)=rprimdtry(:,:) 1322 ! Generate the symmetry operations in the conventional vector coordinates 1323 rprimdconv(:,1)=bravais(3:5) 1324 rprimdconv(:,2)=bravais(6:8) 1325 rprimdconv(:,3)=bravais(9:11) 1326 axes(:,:)=zero 1327 axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one 1328 symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym) 1329 call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym) 1330 if(bravais(1)/=6)then 1331 iaxis=naxes_ortho+1 1332 else 1333 iaxis=naxes_hexa+1 1334 end if 1335 next_stage=0 1336 !DEBUG 1337 ! write(std_out,*)' symbrav: next stage, bravais(1), bravais(2) and symrelconv' 1338 ! write(std_out,'(a,2i4)')' bravais(1:2)=',bravais(1:2) 1339 ! write(std_out,'(a,9f12.6)')' rprimdconv=',rprimdconv(:,:) 1340 ! do isym=1,nsym 1341 ! write(std_out,'(9i4)')symrelconv(:,:,isym) 1342 ! enddo 1343 ! call flush(std_out) 1344 !ENDDEBUG 1345 end if 1346 1347 ! Go to the next iaxis that will be left invariant 1348 iaxis=iaxis-1 1349 do jaxis=iaxis,1,-1 1350 if(bravais(1)/=6)then 1351 axis_trial(:)=ortho_axes(:,jaxis) 1352 else 1353 axis_trial(:)=hexa_axes(:,jaxis) 1354 end if 1355 ! DEBUG 1356 ! write(std_out,*)' symbrav : ixaxis, trial jaxis=',iaxis,jaxis 1357 ! write(std_out,*)' axis_trial=',axis_trial 1358 ! ENDDEBUG 1359 invariant=1 1360 ! Examine whether all symmetry operations leave the axis invariant (might be reversed, though) 1361 do isym=1,nsym 1362 if(sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+(-axis_trial(:))))/=0 .and. & 1363 sum(abs(matmul(symrelconv(:,:,isym),axis_trial)+axis_trial(:)))/=0 )invariant=0 1364 end do 1365 if(invariant==1)then 1366 iaxis=jaxis 1367 !DEBUG 1368 ! write(msg, '(2a,i3)' )ch10,' symbrav : found invariant axis, jaxis=',jaxis 1369 ! call wrtout(std_out,msg) 1370 !ENDDEBUG 1371 exit 1372 end if 1373 end do 1374 1375 if(invariant==0)then 1376 ! Not a single axis was invariant with respect to all operations ?! 1377 ! do isym=1,nsym; write(std_out, '(a,10i4)' )' isym,symrelconv=',isym,symrelconv(:,:,isym); enddo 1378 write(msg, '(3a,3i3,2a,i3,2a,i3)' )& 1379 'Could not succeed to determine the bravais lattice (not a single invariant)',ch10,& 1380 'problem,iaxis,invariant=',problem,iaxis,invariant,ch10,& 1381 'bravais(1)=',bravais(1),ch10,& 1382 'iholohedry=',iholohedry 1383 ABI_BUG(msg) 1384 end if 1385 1386 call matr3inv(rprimdconv,rprimdconv_invt) 1387 axis_red(:)=axis_trial(1)*rprimdconv_invt(1,:)+ & 1388 & axis_trial(2)*rprimdconv_invt(2,:)+ & 1389 & axis_trial(3)*rprimdconv_invt(3,:) 1390 axis_cart(:)=axis_red(1)*rprimdnow(:,1)+ & 1391 & axis_red(2)*rprimdnow(:,2)+ & 1392 & axis_red(3)*rprimdnow(:,3) 1393 norm=sum(axis_cart(:)**2) 1394 !DEBUG 1395 ! write(6,*)' axis_trial =',axis_trial 1396 ! write(6,*)' axis_red =',axis_red 1397 ! write(6,*)' axis_cart =',axis_cart 1398 ! write(6,*)' rprimdnow=',rprimdnow 1399 !ENDDEBUG 1400 ! Expand by a uniform, quite arbitrary, dilatation, along the invariant axis 1401 ! Note : make these dilatation different, according to ideform 1402 ! XG 20151221 : Still, the interplay between the size of the deformation and the tolsym is not easy to address. 1403 ! Indeed the deformation must be sufficiently large to be perceived by symlatt as a real breaking of the 1404 ! symmetry of the lattice. In order to deal with all the small values od tolsym, it has been set at a minimum of tol3, 1405 ! but it must also be larger than tolsym. Moreover, for some axis choice, the deformation is not aligned with the axis, decreasing 1406 ! the effective deformation length. An additional factor of three is thus included, actually increased to six just to be sure... 1407 do ii=1,3 1408 scprod=axis_cart(1)*rprimdnow(1,ii)+axis_cart(2)*rprimdnow(2,ii)+axis_cart(3)*rprimdnow(3,ii) 1409 rprimdtry(:,ii)=rprimdnow(:,ii)+ideform*(max(tol3,six*tolsym)-tol6)*scprod/norm*axis_cart(:) 1410 end do 1411 1412 !DEBUG 1413 ! write(6,*)' rprimdtry=',rprimdtry 1414 !ENDDEBUG 1415 1416 end do ! ideform 1417 1418 if(bravais(1)/=iholohedry)then 1419 write(msg, '(3a,3i3,2a,i3,2a,i3)' )& 1420 'Despite efforts, Could not succeed to determine the bravais lattice :',ch10,& 1421 'bravais(1)=',bravais(1),ch10,& 1422 'iholohedry=',iholohedry 1423 ABI_BUG(msg) 1424 end if 1425 1426 ABI_FREE(symrelconv) 1427 1428 if (PRESENT(axis)) then ! Return symmetry axis. 1429 axis=(/0,0,0/) 1430 if (iaxis/=0) then 1431 if(bravais(1)/=6)then 1432 axis=ortho_axes(:,iaxis) 1433 else 1434 axis=hexa_axes(:,iaxis) 1435 end if 1436 end if 1437 end if 1438 1439 !DEBUG 1440 !write(std_out,'(a)')' symbrav : exit ' 1441 !ENDDEBUG 1442 1443 end subroutine symbrav
m_symfind/symfind [ Functions ]
[ Top ] [ m_symfind ] [ Functions ]
NAME
symfind
FUNCTION
Symmetry finder. From the symmetries of the Bravais lattice (ptsymrel), select those that leave invariant the system, and generate the corresponding tnons vectors and symafm information. Unlike symfind_expert, does NOT resymmetrize atomic positions and tnons for more robust determination of the symmetries. The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972) [[cite:Worlton1972]]
INPUTS
chrgat(natom) (optional)=target charge for each atom. Not always used, it depends on the value of constraint_kind
invardir_red (optional)=reduced coordinates of an invariant direction (only acting with symrel - not tnons)
invar_z (optional)= if 1, the z direction must stay invariant for all symrel applied ;
if 2, z must stay invariant and also there cannot be any associated tnons along z.
gprimd(3,3)=dimensional primitive translations for reciprocal space
msym=default maximal number of symmetries
natom=number of atoms in cell.
nptsym=number of point symmetries of the Bravais lattice
nspden= number of spin-density components. When 4, the three components of spinat are taken into account, instead of only z-component.
nucdipmom(3,natom) (optional) array of nuclear dipole moments
ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
of the Bravais lattice in real space in terms
of primitive translations.
spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
tolsym=tolerance for the symmetries
typat(natom)=integer identifying type of atom.
use_inversion=1 if inversion and improper rotations can be included in set of symmetries
xred(3,natom)=reduced coordinates of atoms in terms of real space
primitive translations
OUTPUT
ierr (optional)=if non-zero, the symmetry operations do not form a group nsym=actual number of symmetries symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations symrel(3,3,1:msym)= nsym symmetry operations in real space in terms of primitive translations tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would be 0 0 0 each for a symmorphic space group)
SOURCE
95 subroutine symfind(gprimd,msym,natom,nptsym,nspden,nsym,& 96 & prtvol, ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,& 97 & chrgat,ierr,nucdipmom,invardir_red,invar_z) ! Optional 98 99 !Arguments ------------------------------------ 100 !scalars 101 integer,intent(in) :: msym,natom,nptsym,nspden,use_inversion 102 integer,intent(in) :: prtvol 103 integer,optional,intent(in) :: invar_z 104 integer,optional,intent(out) :: ierr 105 integer,intent(out) :: nsym 106 real(dp),intent(in) :: tolsym 107 !arrays 108 integer,intent(in) :: ptsymrel(3,3,msym),typat(natom) 109 integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i 110 real(dp),intent(in) :: gprimd(3,3),spinat(3,natom),xred(3,natom) 111 real(dp),optional,intent(in) :: invardir_red(3),chrgat(natom) 112 real(dp),optional, intent(in) :: nucdipmom(3,natom) 113 real(dp),intent(inout) :: tnons(3,msym) !vz_i 114 115 !Local variables------------------------------- 116 !scalars 117 integer :: found3,foundcl,iatom,iatom0,iatom1,iatom2,iatom3,iclass,iclass0,ierr_,ii 118 integer :: isym,jj,kk,natom0,nclass,ntrial,printed,trialafm,trialok 119 real(dp) :: det,diff1,diff2,diff3,diffr1,diffr2,diffr3,ndnorm,nucdipmomcl2,nucdipmomcl20 120 real(dp) :: spinat2,spinatcl2,spinatcl20,tolsym2 121 ! TRUE if antiferro symmetries are used with non-collinear magnetism. 122 integer :: afm_noncoll=1 !For nspden==4. If 1, all symops are permitted ; if 0 symafm must be 1. 123 !For nspden=4. If noncoll_orthorhombic1, require the symmetry operations to be a subset of the orthorhombic symmetries, except if all spinat=0.. 124 integer :: noncoll_orthorhombic=0 125 logical :: test_sameabsspin,test_samechrg 126 logical :: test_samenucdipmom 127 character(len=500) :: msg 128 !arrays 129 integer,allocatable :: class(:,:),natomcl(:),typecl(:) 130 real(dp) :: diff(3),invardir_red_rot(3),hand2(3),hand3(3),ndtest(3),rprimd(3,3),spinat0(3),xred0(3) 131 !real(dp) :: symnucdipmom2(3) 132 real(dp) :: symnucdipmom2cart(3,3),symnucdipmom2red(3,3) 133 real(dp) :: symspinat1(3),symspinat2(3),symxred2(3),trialnons(3) 134 real(dp),allocatable :: chrgat_(:) 135 real(dp),allocatable :: chrgatcl(:) 136 real(dp),allocatable :: local_nucdipmom(:,:,:),nucdipmomcl(:,:),nucdipmomred(:,:,:) 137 real(dp),allocatable :: spinatcl(:,:),spinatred(:,:) 138 139 !************************************************************************** 140 141 !DEBUG 142 !write(std_out,'(a)')' m_symfind%symfind : enter ' 143 !call flush(std_out) 144 !ENDDEBUG 145 146 !DEBUG 147 ! if (prtvol>1) msg="remove me later" 148 ! write(std_out,*)' symfind : enter' 149 ! call flush(6) 150 ! write(std_out,*)' ptsymrel matrices are :' 151 ! do isym=1,nptsym 152 ! write(std_out,'(i4,4x,9i4)' )isym,ptsymrel(:,:,isym) 153 ! end do 154 ! write(std_out,*)' symfind : natom=',natom 155 ! do iatom=1,natom 156 ! write(std_out,*)' atom number',iatom 157 ! write(std_out,*)' typat =',typat(iatom) 158 ! write(std_out,*)' spinat =',spinat(:,iatom) 159 ! write(std_out,*)' xred =',xred(:,iatom) 160 ! end do 161 ! write(std_out,*)' ' 162 ! call flush(6) 163 !ENDDEBUG 164 165 ABI_UNUSED(prtvol) 166 167 !Find the number of classes of atoms (type, chrg and spinat must be identical, 168 !spinat might differ by a sign, if aligned with the z direction, or, 169 ! type and nucdipmom must be identical) 170 !natomcl(iclass) will contain the number of atoms in the class 171 !typecl(iclass) will contain the type of the atoms in the class 172 !chrgcl(iclass) will contain the charge of the atoms in the class 173 !spinatcl(1:3,iclass) will contain the spinat of the atoms in the class 174 !class(1:natomclass(iclass),iclass) will contain the index of the 175 !atoms belonging to the class 176 ABI_MALLOC(class,(natom+3,natom)) 177 ABI_MALLOC(natomcl,(natom)) 178 ABI_MALLOC(typecl,(natom)) 179 ABI_MALLOC(chrgat_,(natom)) 180 ABI_MALLOC(chrgatcl,(natom)) 181 ABI_MALLOC(spinatcl,(3,natom)) 182 ABI_MALLOC(local_nucdipmom,(3,3,natom)) 183 ABI_MALLOC(nucdipmomcl,(3,natom)) 184 185 tolsym2=tolsym**2 186 187 chrgat_(:)=zero 188 if(present(chrgat))then 189 chrgat_(:)=chrgat(:) 190 endif 191 192 local_nucdipmom(:,:,:) = zero 193 if(present(nucdipmom)) then 194 local_nucdipmom(1:3,1,:) = nucdipmom(1:3,:) 195 end if 196 ! for each nuclear dipole we need a local right handed coord system, so we can 197 ! test later for whether a symmetry operation preserves the circulation induced 198 ! by the dipole 199 do iatom=1, natom 200 ndnorm=sqrt(DOT_PRODUCT(local_nucdipmom(1:3,1,iatom),local_nucdipmom(1:3,1,iatom))) 201 202 ! if nuclear dipole has effectively zero size, move on to the next atom 203 if (ndnorm < tol8) cycle 204 205 ! for testing purposes, we care only about direction so renormalize to unity 206 local_nucdipmom(1:3,1,iatom) = local_nucdipmom(1:3,1,iatom)/ndnorm 207 208 ! make a random vector, each component is (0,1] 209 call random_number(ndtest) 210 211 ! vector 2 is constructed to be orthogonal to original nuclear dipole moment vector 212 call acrossb(local_nucdipmom(1:3,1,iatom),ndtest(1:3),local_nucdipmom(1:3,2,iatom)) 213 214 ! vector 3 is orthogonal to 1 and 2, and 1,2,3 form a right-handed set 215 call acrossb(local_nucdipmom(1:3,1,iatom),local_nucdipmom(1:3,2,iatom),local_nucdipmom(1:3,3,iatom)) 216 end do 217 218 ! need rprimd later to transform back to cart coords 219 call matr3inv(gprimd,rprimd) 220 221 !DEBUG 222 !write(std_out,'(a)')' m_symfind%symfind : before initialise with the first atom ' 223 !call flush(std_out) 224 !ENDDEBUG 225 226 !Initialise with the first atom 227 nclass=1 228 natomcl(1)=1 229 typecl(1)=typat(1) 230 chrgatcl(1)=chrgat_(1) 231 spinatcl(:,1)=spinat(:,1) 232 nucdipmomcl(:,1)=local_nucdipmom(:,1,1) 233 class(1,1)=1 234 if(natom>1)then 235 do iatom=2,natom 236 ! DEBUG 237 ! write(std_out,*)' ' 238 ! write(std_out,*)' symfind : examine iatom=',iatom 239 ! ENDDEBUG 240 foundcl=0 241 do iclass=1,nclass 242 ! Compare the typat, chrg and spinat of atom iatom with existing ones. 243 ! At this stage, admit either identical spinat, or spin-flip spinat. 244 if( typat(iatom)==typecl(iclass)) then 245 test_samechrg= (abs(chrgat_(iatom)-chrgatcl(iclass))<tolsym) 246 if(nspden/=4)then 247 test_sameabsspin=(abs(abs(spinat(3,iatom))-abs(spinatcl(3,iclass)))<tolsym) 248 else if(nspden==4)then 249 spinat2 =spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2 250 spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2 251 test_sameabsspin=abs(spinat2-spinatcl2)<tolsym 252 endif 253 test_samenucdipmom= & 254 & abs(local_nucdipmom(1,1,iatom)-nucdipmomcl(1,iclass))<tolsym .and. & 255 & abs(local_nucdipmom(2,1,iatom)-nucdipmomcl(2,iclass))<tolsym .and. & 256 & abs(local_nucdipmom(3,1,iatom)-nucdipmomcl(3,iclass))<tolsym 257 ! note in the following test, m_chkinp/chkinp has already prevented nucdipmom to be 258 ! nonzero when spinat is nonzero 259 if( test_samechrg .and. test_sameabsspin .and. test_samenucdipmom ) then 260 ! DEBUG 261 ! write(std_out,*)' symfind : find it belongs to class iclass=',iclass 262 ! write(std_out,*)' symfind : spinat(:,iatom)=',spinat(:,iatom) 263 ! write(std_out,*)' symfind : spinatcl(:,iclass)=',spinatcl(:,iclass) 264 ! write(std_out,*)' symfind : test_sameabsspin=',test_sameabsspin 265 ! write(std_out,*)' ' 266 ! ENDDEBUG 267 natomcl(iclass)=natomcl(iclass)+1 268 class(natomcl(iclass),iclass)=iatom 269 foundcl=1 270 exit 271 end if 272 end if 273 end do 274 ! If no class with these characteristics exist, create one 275 if(foundcl==0)then 276 nclass=nclass+1 277 natomcl(nclass)=1 278 typecl(nclass)=typat(iatom) 279 chrgatcl(nclass)=chrgat_(iatom) 280 spinatcl(:,nclass)=spinat(:,iatom) 281 nucdipmomcl(:,nclass)=local_nucdipmom(:,1,iatom) 282 class(1,nclass)=iatom 283 end if 284 end do 285 end if 286 287 !DEBUG 288 !write(std_out,*)' ' 289 !write(std_out,*)' symfind : found ',nclass,' nclass of atoms' 290 !do iclass=1,nclass 291 !write(std_out,*)' class number',iclass 292 !write(std_out,*)' natomcl =',natomcl(iclass) 293 !write(std_out,*)' typecl =',typecl(iclass) 294 !write(std_out,*)' spinatcl=',spinatcl(:,iclass) 295 !write(std_out,*)' class =',(class(iatom,iclass),iatom=1,natomcl(iclass)) 296 !end do 297 !write(std_out,*)' ' 298 !ENDDEBUG 299 300 !DEBUG 301 !write(std_out,'(a)')' m_symfind%symfind : before select the class ' 302 !call flush(std_out) 303 !ENDDEBUG 304 305 !Select the class with the least number of atoms, and non-zero spinat if any 306 !It is important to select a magnetic class of atom, if any, otherwise 307 !the determination of the initial (inclusive) set of symmetries takes only 308 !non-magnetic symmetries, and not both magnetic and non-magnetic ones, see later. 309 !On the contrary, the chrgat_ data does not play any role, it is invariant upon atomic-centered symmetries 310 iclass0=1 311 natom0=natomcl(1) 312 spinatcl20=spinatcl(1,1)**2+spinatcl(2,1)**2+spinatcl(3,1)**2 313 nucdipmomcl20=nucdipmomcl(1,1)**2+nucdipmomcl(2,1)**2+nucdipmomcl(3,1)**2 314 if(nclass>1)then 315 do iclass=2,nclass 316 spinatcl2=spinatcl(1,iclass)**2+spinatcl(2,iclass)**2+spinatcl(3,iclass)**2 317 nucdipmomcl2=nucdipmomcl(1,iclass)**2+nucdipmomcl(2,iclass)**2+nucdipmomcl(3,iclass)**2 318 if( (natomcl(iclass)<natom0 & 319 & .and. .not. (spinatcl20>tolsym .and. spinatcl2<tolsym) & 320 & .and. .not. (nucdipmomcl20>tolsym .and. nucdipmomcl2<tolsym) ) & 321 & .or. (spinatcl20<tolsym .and. spinatcl2>tolsym) & 322 & .or. (nucdipmomcl20<tolsym .and. nucdipmomcl2>tolsym) )then 323 iclass0=iclass 324 natom0=natomcl(iclass) 325 spinatcl20=spinatcl2 326 nucdipmomcl20=nucdipmomcl2 327 end if 328 end do 329 end if 330 331 printed=0 332 333 !If non-collinear spinat have to be used, transfer them in reduced coordinates 334 if (nspden==4) then 335 ABI_MALLOC(spinatred,(3,natom)) 336 do iatom=1,natom 337 do ii=1,3 338 spinatred(1:3,iatom)=MATMUL(TRANSPOSE(gprimd),spinat(1:3,iatom)) 339 end do 340 end do 341 end if 342 343 !DEBUG 344 !write(std_out,*)' ' 345 !write(std_out,*)' symfind : has selected iclass0=',iclass0 346 !write(std_out,*)' # iatom xred spinat (spinatred if nspden=4) ' 347 !do iatom0=1,natomcl(iclass0) 348 !iatom=class(iatom0,iclass0) 349 !if(nspden/=4)then 350 ! write(std_out,'(2i4,6f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom) 351 !else if(nspden==4)then 352 ! write(std_out,'(2i4,9f10.4)' )iatom0,iatom,xred(:,iatom),spinat(:,iatom),spinatred(:,iatom) 353 !endif 354 !end do 355 !write(std_out,*)' ' 356 !ENDDEBUG 357 358 359 !represent nuclear dipole moments in reduced coords 360 ABI_MALLOC(nucdipmomred,(3,3,natom)) 361 do iatom=1,natom 362 do ii=1,3 363 nucdipmomred(1:3,ii,iatom)=MATMUL(TRANSPOSE(gprimd),local_nucdipmom(1:3,ii,iatom)) 364 end do 365 end do 366 367 !DEBUG 368 !write(std_out,'(a)')' m_symfind%symfind : before big loop ' 369 !call flush(std_out) 370 !ENDDEBUG 371 372 !Big loop over each symmetry operation of the Bravais lattice 373 nsym=0 374 do isym=1,nptsym 375 376 !DEBUG 377 !write(std_out,'(a,i4)')' m_symfind%symfind : enter loop isym=',isym 378 !call flush(std_out) 379 !ENDDEBUG 380 381 if(present(invardir_red))then 382 ! ji: Check whether symmetry operation leaves invardir_red invariant 383 invardir_red_rot(:) = ptsymrel(:,1,isym)*invardir_red(1) + & 384 & ptsymrel(:,2,isym)*invardir_red(2) + & 385 & ptsymrel(:,3,isym)*invardir_red(3) 386 diff(:)=invardir_red(:)-invardir_red_rot(:) 387 if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 ) cycle 388 endif 389 390 !DEBUG 391 !write(std_out,'(a,i4)')' m_symfind%symfind : 1' 392 !call flush(std_out) 393 !ENDDEBUG 394 395 if (use_inversion==0) then 396 det=ptsymrel(1,1,isym)*ptsymrel(2,2,isym)*ptsymrel(3,3,isym)+& 397 & ptsymrel(2,1,isym)*ptsymrel(3,2,isym)*ptsymrel(1,3,isym)+& 398 & ptsymrel(1,2,isym)*ptsymrel(2,3,isym)*ptsymrel(3,1,isym) - & 399 & (ptsymrel(3,1,isym)*ptsymrel(2,2,isym)*ptsymrel(1,3,isym)+& 400 & ptsymrel(2,1,isym)*ptsymrel(1,2,isym)*ptsymrel(3,3,isym)+& 401 & ptsymrel(3,2,isym)*ptsymrel(2,3,isym)*ptsymrel(1,1,isym)) 402 if(det==-1) cycle 403 end if 404 405 !DEBUG 406 !write(std_out,'(a,i4)')' m_symfind%symfind : 2' 407 !call flush(std_out) 408 !ENDDEBUG 409 410 ! jellium slab and spatially varying chemical potential cases: 411 ! (actually, an inversion symmetry/mirror plane perpendicular to z symmetry operation might still be allowed... TO BE DONE !) 412 if(present(invar_z))then 413 if (invar_z/=0) then 414 ! check whether symmetry operation produce a rotation only in the xy plane 415 if( ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. & 416 & ptsymrel(3,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle 417 ! check whether symmetry operation does not change the z 418 if( ptsymrel(3,3,isym)/=1 ) cycle 419 end if 420 end if 421 422 !DEBUG 423 !write(std_out,'(a,i4)')' m_symfind%symfind : 3' 424 !call flush(std_out) 425 !ENDDEBUG 426 427 ! If noncoll_orthorhombic=1, require orthorhombic operations of symmetries, except if spinat=0. 428 if (nspden==4 .and. noncoll_orthorhombic==1)then 429 if(sum(abs(spinat(:,:)))>tol14)then 430 if( ptsymrel(1,3,isym)/=0 .or. ptsymrel(2,3,isym)/=0 .or. & 431 & ptsymrel(1,2,isym)/=0 .or. ptsymrel(3,2,isym)/=0 .or. & 432 & ptsymrel(2,1,isym)/=0 .or. ptsymrel(3,2,isym)/=0 ) cycle 433 endif 434 endif 435 436 !DEBUG 437 !write(std_out,'(a,i4)')' m_symfind%symfind : 4' 438 !call flush(std_out) 439 !ENDDEBUG 440 441 ! Select a tentative set of associated translations 442 ! First compute the symmetric of the first atom in the smallest class, 443 ! using the point symmetry, and also the symmetric of spinat(red). 444 iatom0=class(1,iclass0) 445 xred0(:)=ptsymrel(:,1,isym)*xred(1,iatom0)+ & 446 & ptsymrel(:,2,isym)*xred(2,iatom0)+ & 447 & ptsymrel(:,3,isym)*xred(3,iatom0) 448 if (nspden/=4) then 449 spinat0(:)=spinat(:,iatom0) 450 else 451 spinat0(:)=ptsymrel(:,1,isym)*spinatred(1,iatom0)+ & 452 & ptsymrel(:,2,isym)*spinatred(2,iatom0)+ & 453 & ptsymrel(:,3,isym)*spinatred(3,iatom0) 454 endif 455 456 !DEBUG 457 !write(std_out,'(a,i4)')' m_symfind%symfind : 5' 458 !call flush(std_out) 459 !ENDDEBUG 460 461 ! From the set of possible images, deduce tentative translations, 462 ! and magnetic factor then test whether it send each atom on a symmetric one 463 ntrial=0 464 do ii=1,natom0 465 !DEBUG 466 ! write(std_out,'(a,2i4)')' symfind : loop isym,ii=',isym,ii 467 !ENDDEBUG 468 iatom1=class(ii,iclass0) 469 470 ! The tentative translation is found 471 trialnons(:)=xred(:,iatom1)-xred0(:) 472 ! Compare the spinat vectors 473 if (nspden/=4) then 474 symspinat1(:)=spinat(:,iatom1) 475 else 476 symspinat1(:)=spinatred(:,iatom1) 477 end if 478 479 !DEBUG 480 ! write(std_out,'(a,6f10.4)')' symspinat1,spinat0=',symspinat1(:),spinat0(:) 481 !ENDDEBUG 482 483 trialafm=1 484 if(sum(abs(symspinat1(:)-spinat0(:)))>tolsym)then 485 trialafm=-1 486 if(nspden==4 .and. afm_noncoll==0)cycle 487 if(sum(abs(symspinat1(:)+spinat0(:)))>tolsym)cycle 488 endif 489 490 if(sum(abs(local_nucdipmom(:,1,iatom1)-local_nucdipmom(:,1,iatom0)))>tolsym)then 491 write(msg,'(3a,3i5)')& 492 'Problem with matching the nuclear dipole moment within a class.',ch10,& 493 'isym,iatom0,iatom1=',isym,iatom0,iatom1 494 ABI_ERROR(msg) 495 end if 496 497 ! jellium slab case: check whether symmetry operation has no translational 498 ! component along z 499 if(present(invar_z))then 500 if( invar_z==2 .and. abs(trialnons(3)) > tolsym ) cycle 501 endif 502 trialok=1 503 504 ! DEBUG 505 ! write(std_out, '(a,i3,a,i3,a,i3,a,3f12.4,i3)') ' Try isym=',isym,' sending iatom0 ',iatom0,' to iatom1 ',iatom1,' with trialnons(:),trialafm =',trialnons(:),trialafm 506 ! ENDDEBUG 507 508 ! Loop over all classes, then all atoms in the class, 509 ! to find whether they have a symmetric 510 do iclass=1,nclass 511 do jj=1,natomcl(iclass) 512 513 iatom2=class(jj,iclass) 514 ! Generate the tentative symmetric position of iatom2 515 symxred2(:)=ptsymrel(:,1,isym)*xred(1,iatom2)+ & 516 & ptsymrel(:,2,isym)*xred(2,iatom2)+ & 517 & ptsymrel(:,3,isym)*xred(3,iatom2)+ trialnons(:) 518 ! Generate the tentative symmetric spinat of iatom2 519 if (nspden/=4) then 520 symspinat2(:)=trialafm*spinat(:,iatom2) 521 else 522 symspinat2(:)=trialafm*(ptsymrel(:,1,isym)*spinatred(1,iatom2)+ & 523 & ptsymrel(:,2,isym)*spinatred(2,iatom2)+ & 524 & ptsymrel(:,3,isym)*spinatred(3,iatom2)) 525 end if 526 ! Generate the tentative symmetric nucdipmom of iatom2 527 do kk = 1, 3 528 symnucdipmom2red(:,kk)=ptsymrel(:,1,isym)*nucdipmomred(1,kk,iatom2)+ & 529 & ptsymrel(:,2,isym)*nucdipmomred(2,kk,iatom2)+ & 530 & ptsymrel(:,3,isym)*nucdipmomred(3,kk,iatom2) 531 ! transform back to cart coords for final comparison to nucdipmom 532 symnucdipmom2cart(:,kk)=MATMUL(rprimd,symnucdipmom2red(:,kk)) 533 end do 534 535 ! DEBUG 536 ! write(std_out,'(a,i4,a,3f8.4,a,3f8.4,a,3f8.4)')& 537 !& ' Test iatom2=',iatom2,' at xred=',xred(:,iatom2),'. Is sent to',symxred2(:),' with symspinat2=',symspinat2(:) 538 ! ENDDEBUG 539 540 ! Check whether there exists an atom of the same class at the 541 ! same location, with the correct spinat and nuclear dipole moment circulation 542 do kk=1,natomcl(iclass) 543 544 found3=1 545 iatom3=class(kk,iclass) 546 ! Check the location 547 diffr1=xred(1,iatom3)-symxred2(1) 548 diff1=diffr1-nint(diffr1) 549 if(diff1**2>tolsym2)then 550 found3=0 ; cycle 551 else 552 diffr2=xred(2,iatom3)-symxred2(2) 553 diff2=diffr2-nint(diffr2) 554 if(diff2**2>tolsym2)then 555 found3=0 ; cycle 556 else 557 diffr3=xred(3,iatom3)-symxred2(3) 558 diff3=diffr3-nint(diffr3) 559 if( (diff1**2+diff2**2+diff3**2) > tolsym**2 )then 560 found3=0 ; cycle 561 endif 562 endif 563 endif 564 ! Check the spinat 565 if (nspden/=4) then 566 diff(:)=spinat(:,iatom3)-symspinat2(:) 567 else 568 diff(:)=spinatred(:,iatom3)-symspinat2(:) 569 end if 570 if( (diff(1)**2+diff(2)**2+diff(3)**2) > tolsym**2 )then 571 found3=0 572 cycle 573 endif 574 ! Check the nucdipmom 575 ! hand3 gives original circulation sense of nuclear dipole 576 call acrossb(local_nucdipmom(1:3,2,iatom3),local_nucdipmom(1:3,3,iatom3),hand3) 577 578 ! hand2 gives circulation sense of tentative, symmetry equivalent nuclear dipole 579 call acrossb(symnucdipmom2cart(1:3,2),symnucdipmom2cart(1:3,3),hand2) 580 581 diff(:)=hand3(:)-hand2(:) 582 if( any(abs(diff)>tolsym) )found3=0 583 584 if(found3==1)exit 585 end do ! End loop over iatom3 586 587 if(found3==0)then 588 trialok=0 589 exit 590 end if 591 end do ! End loop over iatom2 592 593 if(trialok==0)exit 594 end do ! End loop over all classes 595 596 !DEBUG 597 ! write(std_out,*)' For trial isym=',isym,', trialok = ',trialok 598 ! write(std_out,*)' ' 599 !ENDDEBUG 600 601 if(trialok==1)then 602 nsym=nsym+1 603 if(nsym>msym)then 604 write(msg,'(a,i0,2a,i0,4a)')& 605 'The number of symmetries (including non-symmorphic translations) is: ', nsym, ch10,& 606 'is larger than maxnsym: ',msym,ch10,& 607 'Action: increase maxnsym in the input, or take a cell that is primitive, ',ch10,& 608 'or at least smaller than the present one.' 609 ABI_ERROR(msg) 610 end if 611 ntrial=ntrial+1 612 symrel(:,:,nsym)=ptsymrel(:,:,isym) 613 symafm(nsym)=trialafm 614 tnons(:,nsym)=trialnons(:)-nint(trialnons(:)-tolsym) 615 end if 616 617 end do ! End the loop on tentative translations 618 end do ! End big loop over each symmetry operation of the Bravais lattice 619 620 !DEBUG 621 !write(std_out,'(a)')' m_symfind%symfind : after big loop, will call ABI_FREE ' 622 !call flush(std_out) 623 !ENDDEBUG 624 625 ABI_FREE(class) 626 ABI_FREE(natomcl) 627 ABI_FREE(chrgat_) 628 ABI_FREE(chrgatcl) 629 ABI_FREE(spinatcl) 630 ABI_FREE(typecl) 631 ABI_FREE(local_nucdipmom) 632 ABI_FREE(nucdipmomcl) 633 ABI_FREE(nucdipmomred) 634 if (nspden==4) then 635 ABI_FREE(spinatred) 636 end if 637 638 !DEBUG 639 !write(std_out,'(a,i6)')' m_symfind%symfind : call sg_multable, nsym= ',nsym 640 !call flush(std_out) 641 !ENDDEBUG 642 643 ! The algorithm in sg_multable is still cubic in nsym, so avoid calling it uselessly when nsym is too large 644 if(present(ierr) .or. nsym<=384)then 645 call sg_multable(nsym, symafm, symrel, ierr_, tnons=tnons, tnons_tol=tolsym) 646 else 647 ierr_=0 648 endif 649 650 !DEBUG 651 !write(std_out,'(a)')' m_symfind%symfind : call print_symmetries, ierr_= ',ierr_ 652 !call flush(std_out) 653 !ENDDEBUG 654 655 if (ierr_/=0) then 656 call print_symmetries(nsym,symrel,tnons,symafm) 657 end if 658 659 if(.not.present(ierr))then 660 ABI_CHECK(ierr_==0,"Error in group closure") 661 else 662 ierr=ierr_ 663 endif 664 665 !DEBUG 666 ! write(msg,'(a,I0,es16.6,a)')' symfind : exit, nsym, tolsym=',nsym,tolsym,ch10 667 ! write(msg,'(2a)') trim(msg),' symrel matrices, symafm and tnons are :' 668 ! call wrtout(std_out,msg) 669 ! do isym=1,nsym 670 ! write(msg,'(i4,4x,3i4,2x,3i4,2x,3i4,4x,i4,4x,3f8.4)' ) isym,symrel(:,:,isym),& 671 ! & symafm(isym),tnons(:,isym) 672 ! call wrtout(std_out,msg) 673 ! end do 674 !stop 675 !ENDDEBUG 676 677 !DEBUG 678 !write(std_out,'(a)')' m_symfind%symfind : exit ' 679 !call flush(std_out) 680 !ENDDEBUG 681 682 end subroutine symfind
m_symfind/symfind_expert [ Functions ]
[ Top ] [ m_symfind ] [ Functions ]
NAME
symfind_expert
FUNCTION
Symmetry finder, with an added layer of robustness compared to symfind, and for which resymmetrization of atomic positions and tnons is needed.. From the symmetries of the Bravais lattice (ptsymrel), select those that leave invariant the system, and generate the corresponding tnons vectors and symafm information. Unlike symfind_expert, does NOT resymmetrize atomic positions and tnons for more robust determination of the symmetries. The algorithm is explained in T.G. Worlton and J.L. Warren, Comp. Phys. Comm. 3, 88 (1972) [[cite:Worton1972]]
INPUTS
chrgat(natom) (optional)=target charge for each atom. Not always used, it depends on the value of constraint_kind
invardir_red (optional)=reduced coordinates of an invariant direction (only acting with symrel - not tnons)
invar_z (optional)= if 1, the z direction must stay invariant for all symrel applied ;
if 2, z must stay invariant and also there cannot be any associated tnons along z.
gprimd(3,3)=dimensional primitive translations for reciprocal space
msym=default maximal number of symmetries
natom=number of atoms in cell.
nptsym=number of point symmetries of the Bravais lattice
nspden= number of spin-density components. When 4, the three components of spinat are taken into account, instead of only z-component.
nucdipmom(3,natom) (optional) array of nuclear dipole moments
ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
of the Bravais lattice in real space in terms
of primitive translations.
spinat(3,natom)=initial spin of each atom, in unit of hbar/2.
tolsym=tolerance for the symmetries
typat(natom)=integer identifying type of atom.
usepaw= 0 for non paw calculation; =1 for paw calculation
OUTPUT
nsym=actual number of symmetries symafm(1:msym)=(anti)ferromagnetic part of nsym symmetry operations symrel(3,3,1:msym)= nsym symmetry operations in real space in terms of primitive translations tnons(3,1:msym)=nonsymmorphic translations for each symmetry (would be 0 0 0 each for a symmorphic space group)
SIDE EFFECTS
xred(3,natom)=reduced coordinates of atoms in terms of real space primitive translations. Might be changed during the resymmetrization.
SOURCE
731 subroutine symfind_expert(gprimd,msym,natom,nptsym,nspden,nsym,& 732 & pawspnorb,prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,usepaw,xred,& 733 & chrgat,nucdipmom,invardir_red,invar_z) ! Optional - although for the time being all are required ... 734 735 !Arguments ------------------------------------ 736 !scalars 737 integer,intent(in) :: msym,natom,nptsym,nspden,pawspnorb,usepaw 738 integer,intent(in) :: prtvol 739 integer,optional,intent(in) :: invar_z 740 integer,intent(out) :: nsym 741 real(dp),intent(in) :: tolsym 742 !arrays 743 integer,intent(in) :: ptsymrel(3,3,msym),typat(natom) 744 integer,intent(inout) :: symafm(msym),symrel(3,3,msym) !vz_i 745 real(dp),intent(in) :: gprimd(3,3),spinat(3,natom) 746 real(dp),intent(inout) :: xred(3,natom) 747 real(dp),optional,intent(in) :: invardir_red(3),chrgat(natom) 748 real(dp),optional, intent(in) :: nucdipmom(3,natom) 749 real(dp),intent(inout) :: tnons(3,msym) !vz_i 750 751 !Local variables------------------------------- 752 !scalars 753 integer, save :: print_comment_tolsym=1 754 integer :: fixed_mismatch,mismatch_fft_tnons 755 integer :: ierr,isym,use_inversion 756 character(len=1000) :: msg 757 !arrays 758 integer,allocatable :: indsym(:,:,:),symrec(:,:,:) 759 real(dp),allocatable :: tnons_new(:,:) 760 761 762 !************************************************************************** 763 764 !DEBUG 765 ! write(std_out,*)' m_symfind%symfind_expert : enter ' 766 !ENDDEBUG 767 768 use_inversion=1 769 if (usepaw == 1 .and. (nspden==4.or.pawspnorb>0)) then 770 ABI_COMMENT("Removing inversion and improper rotations from initial space group because of PAW + SOC") 771 use_inversion=0 772 end if 773 774 !DEBUG 775 ! write(std_out,*)' m_symfind%symfind_expert : before call symfind (1) ' 776 !ENDDEBUG 777 778 call symfind(gprimd,msym,natom,nptsym,nspden,nsym,& 779 prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,& 780 chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invar_z=invar_z) 781 782 !DEBUG 783 ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (1) ' 784 !ENDDEBUG 785 786 !If the group closure is not obtained, which should be exceptional, try with a larger tolsym (three times larger) 787 if(ierr/=0)then 788 ABI_WARNING('Will try to obtain group closure by using a tripled tolsym.') 789 call symfind(gprimd,msym,natom,nptsym,nspden,nsym,& 790 prtvol,ptsymrel,spinat,symafm,symrel,tnons,three*tolsym,typat,use_inversion,xred,& 791 chrgat=chrgat,nucdipmom=nucdipmom,ierr=ierr,invardir_red=invardir_red,invar_z=invar_z) 792 ABI_CHECK(ierr==0,"Error in group closure") 793 ABI_WARNING('Succeeded to obtain group closure by using a tripled tolsym.') 794 endif 795 796 ! If the tolerance on symmetries is bigger than 1.e-8, symmetrize tnons for gliding or screw operations, 797 ! symmetrize the atomic positions and recompute the symmetry operations 798 if(tolsym>1.00001e-8)then 799 800 call symmetrize_tnons(nsym,symrel,tnons,tolsym) 801 ABI_MALLOC(indsym,(4,natom,nsym)) 802 ABI_MALLOC(symrec,(3,3,nsym)) 803 do isym=1,nsym 804 call mati3inv(symrel(:,:,isym),symrec(:,:,isym)) 805 end do 806 call symatm(indsym,natom,nsym,symrec,tnons,tolsym,typat,xred) 807 call symmetrize_xred(natom,nsym,symrel,tnons,xred,indsym=indsym) 808 ABI_FREE(indsym) 809 ABI_FREE(symrec) 810 811 if(print_comment_tolsym==1)then 812 write(msg,'(a,es12.3,18a)')& 813 & 'The tolerance on symmetries =',tolsym,' is bigger than 1.0e-8.',ch10,& 814 & 'In order to avoid spurious effects, the atomic coordinates have been',ch10,& 815 & 'symmetrized before storing them in the dataset internal variable.',ch10,& 816 & 'So, do not be surprised by the fact that your input variables (xcart, xred, ...)',ch10,& 817 & 'do not correspond exactly to the ones echoed by ABINIT, the latter being used to do the calculations.',ch10,& 818 & 'This is not a problem per se.',ch10,& 819 & 'Still, in order to avoid this symmetrization (e.g. for specific debugging/development),',& 820 & ' decrease tolsym to 1.0e-8 or lower,',ch10,& 821 & 'or (much preferred) use input primitive vectors that are accurate to better than 1.0e-8.',ch10,& 822 & 'This message will only be printed once, even if there are other datasets where tolsym is bigger than 1.0e-8.' 823 ABI_COMMENT(msg) 824 print_comment_tolsym=0 825 endif 826 827 !DEBUG 828 !write(std_out,*)' m_symfind%symfind_expert : before call symfind (3) ' 829 !ENDDEBUG 830 831 call symfind(gprimd,msym,natom,nptsym,nspden,nsym,& 832 prtvol,ptsymrel,spinat,symafm,symrel,tnons,tolsym,typat,use_inversion,xred,& 833 chrgat=chrgat,nucdipmom=nucdipmom,invardir_red=invardir_red,invar_z=invar_z) 834 835 !DEBUG 836 ! write(std_out,*)' m_symfind%symfind_expert : after call symfind (3) ' 837 !ENDDEBUG 838 839 !Needs one more resymmetrization, for the tnons 840 ABI_MALLOC(tnons_new,(3,nsym)) 841 842 call symmetrize_xred(natom,nsym,symrel,tnons,xred,& 843 & fixed_mismatch=fixed_mismatch,mismatch_fft_tnons=mismatch_fft_tnons,tnons_new=tnons_new,tolsym=tolsym) 844 tnons(:,1:nsym)=tnons_new(:,:) 845 ABI_FREE(tnons_new) 846 end if ! tolsym >1.00001e-8 847 848 !DEBUG 849 ! write(std_out,*)' m_symfind%symfind_expert : exit ' 850 !ENDDEBUG 851 852 end subroutine symfind_expert
m_symfind/symlatt [ Functions ]
[ Top ] [ m_symfind ] [ Functions ]
NAME
symlatt
FUNCTION
From the unit cell vectors (rprimd) and the corresponding metric tensor, find the Bravais lattice and its symmetry operations (ptsymrel). 1) Find the shortest possible primitive vectors for the lattice 2) Determines the holohedral group of the lattice, and the axes to be used for the conventional cell (this is a delicate part, in which the centering of the reduced cell must be taken into account) The idea is to determine the basis vectors of the conventional cell from the reduced cell basis vectors. 3) Generate the symmetry operations of the holohedral group
INPUTS
iout=unit number of output file msym=default maximal number of symmetries. WARNING : cannot be simply set to nsym, because the number of symmetries found here will likely be bigger than sym ! rprimd(3,3)=dimensional primitive translations for real space (bohr) tolsym=tolerance for the symmetries
OUTPUT
bravais(11): bravais(1)=iholohedry
bravais(2)=center
bravais(3:11)=coordinates of rprim in the axes
of the conventional bravais lattice (*2 if center/=0)
nptsym=number of point symmetries of the Bravais lattice
ptsymrel(3,3,1:msym)= nptsym point-symmetry operations
of the Bravais lattice in real space in terms of primitive translations.
NOTES
WARNING: bravais(1) might be given a negative value in another routine, if the cell is non-primitive. The holohedral groups are numbered as follows (see international tables for crystallography (1983), p. 13) iholohedry=1 triclinic 1bar iholohedry=2 monoclinic 2/m iholohedry=3 orthorhombic mmm iholohedry=4 tetragonal 4/mmm iholohedry=5 trigonal 3bar m iholohedry=6 hexagonal 6/mmm iholohedry=7 cubic m3bar m Centering center=0 no centering center=-1 body-centered center=-3 face-centered center=1 A-face centered center=2 B-face centered center=3 C-face centered
SOURCE
1920 subroutine symlatt(bravais,iout,msym,nptsym,ptsymrel,rprimd,tolsym) 1921 1922 !Arguments ------------------------------------ 1923 !scalars 1924 integer,intent(in) :: iout,msym 1925 integer,intent(out) :: nptsym 1926 real(dp),intent(in) :: tolsym 1927 !arrays 1928 integer,intent(out) :: bravais(11),ptsymrel(3,3,msym) 1929 real(dp),intent(in) :: rprimd(3,3) 1930 1931 !Local variables------------------------------- 1932 !scalars 1933 integer,parameter :: mgen=4 1934 integer :: center,fact,found,foundc,ia,iaxis1,iaxis2 1935 integer :: isign1,isign2,ib,icase,igen,iholohedry,ii,index,isym 1936 integer :: itrial,jj,jsym,ngen=0,orthogonal,sign12,sign13,sign23,sumsign 1937 real(dp) :: determinant,norm2a,norm2b,norm2c,norm2trial,reduceda,reducedb,sca 1938 real(dp) :: scalarprod,scb,trace,trace_best,val 1939 character(len=500) :: msg 1940 !arrays 1941 integer,parameter :: list_holo(7)=(/7,6,4,3,5,2,1/) 1942 integer :: ang90(3),equal(3),gen(3,3,mgen),gen2xy(3,3),gen2y(3,3),gen2z(3,3) 1943 integer :: gen3(3,3),gen6(3,3),icoord(3,3),identity(3,3),nvecta(3),nvectb(3) 1944 integer :: order(mgen) 1945 real(dp) :: axes(3,3),axesinvt(3,3),axes_best(3,3),axes_try(3,3) 1946 real(dp) :: cell_base(3,3),coord(3,3),metmin(3,3) 1947 real(dp) :: minim(3,3),scprods(3,3),vecta(3),vectb(3),vectc(3),vin1(3),vin2(3),vext(3) 1948 1949 !************************************************************************** 1950 1951 !DEBUG 1952 !write(std_out,'(a,es14.6)') ' m_symfind%symlatt : enter, tolsym= ',tolsym 1953 !call flush(std_out) 1954 !ENDDEBUG 1955 1956 identity(:,:)=0 ; identity(1,1)=1 ; identity(2,2)=1 ; identity(3,3)=1 1957 nvecta(1)=2 ; nvectb(1)=3 1958 nvecta(2)=1 ; nvectb(2)=3 1959 nvecta(3)=1 ; nvectb(3)=2 1960 1961 !-------------------------------------------------------------------------- 1962 !Reduce the input vectors to a set of minimal vectors 1963 call smallprim(metmin,minim,rprimd) 1964 1965 !DEBUG 1966 !write(std_out,*)' symlatt : minim(:,1)=',minim(:,1) 1967 !write(std_out,*)' symlatt : minim(:,2)=',minim(:,2) 1968 !write(std_out,*)' symlatt : minim(:,3)=',minim(:,3) 1969 !call flush(std_out) 1970 !ENDDEBUG 1971 1972 !-------------------------------------------------------------------------- 1973 !Examine the angles and vector lengths 1974 ang90(:)=0 1975 if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1 1976 if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1 1977 if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1 1978 equal(:)=0 1979 if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1 1980 if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1 1981 if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1 1982 1983 !DEBUG 1984 !write(std_out,*)' ang90=',ang90(:) 1985 !write(std_out,*)' equal=',equal(:) 1986 !call flush(std_out) 1987 !ENDDEBUG 1988 1989 !----------------------------------------------------------------------- 1990 !Identification of the centering 1991 1992 foundc=0 1993 !Default values 1994 fact=1 ; center=0 1995 cell_base(:,:)=minim(:,:) 1996 1997 !Examine each holohedral group 1998 !This search is ordered : should not be happy with tetragonal, 1999 !while there is FCC ... 2000 do index=1,6 2001 2002 ! If the holohedry is already found, exit 2003 if(foundc==1)exit 2004 2005 ! Initialize the target holohedry 2006 iholohedry=list_holo(index) 2007 2008 ! DEBUG 2009 ! write(std_out,*) 2010 ! write(std_out,*)' symlatt : trial holohedry',iholohedry 2011 ! ENDDEBUG 2012 2013 orthogonal=0 2014 if(iholohedry==7 .or. iholohedry==4 .or. iholohedry==3)orthogonal=1 2015 2016 ! Now, will examine different working hypothesis. 2017 ! The set of these hypothesis is thought to cover all possible cases ... 2018 2019 ! Working hypothesis : the basis is orthogonal 2020 if(ang90(1)+ang90(2)+ang90(3)==3 .and. orthogonal==1)then 2021 fact=1 ; center=0 2022 cell_base(:,:)=minim(:,:) 2023 ! Checks that the basis vectors are OK for the target holohedry 2024 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2025 end if 2026 2027 ! Select one trial direction 2028 do itrial=1,3 2029 2030 ! If the holohedry is already found, exit 2031 if(foundc==1)exit 2032 2033 ia=nvecta(itrial) ; ib=nvectb(itrial) 2034 2035 ! This is in case of hexagonal holohedry 2036 if(foundc==0 .and. iholohedry==6 .and. ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then 2037 reduceda=metmin(ib,ia)/metmin(ia,ia) 2038 fact=1 ; center=0 2039 if(abs(reduceda+0.5d0)<tolsym)then 2040 cell_base(:,1)=minim(:,ia) 2041 cell_base(:,2)=minim(:,ib) 2042 cell_base(:,3)=minim(:,itrial) 2043 !DEBUG 2044 ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1) 2045 ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2) 2046 ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3) 2047 !ENDDEBUG 2048 ! Checks that the basis vectors are OK for the target holohedry 2049 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2050 else if(abs(reduceda-0.5d0)<tolsym)then 2051 cell_base(:,1)=minim(:,ia) 2052 cell_base(:,2)=minim(:,ib)-minim(:,ia) 2053 cell_base(:,3)=minim(:,itrial) 2054 ! Checks that the basis vectors are OK for the target holohedry 2055 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2056 end if 2057 end if 2058 2059 ! Working hypothesis : the conventional cell is orthogonal, 2060 ! and the two other vectors are axes of the conventional cell 2061 if(foundc==0 .and. orthogonal==1 .and. ang90(itrial)==1)then 2062 2063 ! Compute the reduced coordinate of trial vector in the basis 2064 ! of the two other vectors 2065 reduceda=metmin(itrial,ia)/metmin(ia,ia) 2066 reducedb=metmin(itrial,ib)/metmin(ib,ib) 2067 cell_base(:,ia)=minim(:,ia) 2068 cell_base(:,ib)=minim(:,ib) 2069 if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. & 2070 & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then 2071 if(abs(abs(reduceda)-0.5d0)<tolsym)center=ib 2072 if(abs(abs(reducedb)-0.5d0)<tolsym)center=ia 2073 fact=2 2074 cell_base(:,itrial)= & 2075 & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0 2076 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2077 else if( abs(abs(reduceda)-0.5d0)<tolsym .and.& 2078 & abs(abs(reducedb)-0.5d0)<tolsym ) then 2079 fact=2 ; center=-1 2080 cell_base(:,itrial)= & 2081 & (minim(:,itrial)-reduceda*minim(:,ia)-reducedb*minim(:,ib) )*2.0d0 2082 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2083 end if 2084 end if 2085 2086 ! Working hypothesis : the conventional cell is orthogonal, and 2087 ! the trial vector is one of the future axes, and the face perpendicular to it is centered 2088 if(foundc==0 .and. iholohedry==3 .and. & 2089 & ang90(ia)==1 .and. ang90(ib)==1 .and. equal(itrial)==1 )then 2090 fact=2 ; center=itrial 2091 cell_base(:,ia)=minim(:,ia)+minim(:,ib) 2092 cell_base(:,ib)=minim(:,ia)-minim(:,ib) 2093 cell_base(:,itrial)=minim(:,itrial) 2094 ! Checks that the basis vectors are OK for the target holohedry 2095 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2096 end if 2097 2098 ! DEBUG 2099 ! write(std_out,*)' after test_b, foundc=',foundc 2100 ! ENDDEBUG 2101 2102 ! Working hypothesis : the conventional cell is orthogonal, and 2103 ! the trial vector is one of the future axes 2104 if(foundc==0 .and. orthogonal==1)then 2105 ! Compute the projection of the two other vectors on the trial vector 2106 reduceda=metmin(itrial,ia)/metmin(itrial,itrial) 2107 reducedb=metmin(itrial,ib)/metmin(itrial,itrial) 2108 ! If both projections are half-integer, one might have found an axis 2109 if( abs(abs(reduceda)-0.5d0)<tolsym .and.& 2110 & abs(abs(reducedb)-0.5d0)<tolsym ) then 2111 vecta(:)=minim(:,ia)-reduceda*minim(:,itrial) 2112 vectb(:)=minim(:,ib)-reducedb*minim(:,itrial) 2113 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2114 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 2115 scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3) 2116 ! Note the order of selection : body-centered is prefered 2117 ! over face centered, which is correct for the tetragonal case 2118 if(abs(norm2a-norm2b)<tolsym*half*(norm2a+norm2b))then 2119 ! The lattice is body centered 2120 fact=2 ; center=-1 2121 cell_base(:,ia)=vecta(:)+vectb(:) 2122 cell_base(:,ib)=vecta(:)-vectb(:) 2123 cell_base(:,itrial)=minim(:,itrial) 2124 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2125 else if(abs(scalarprod)<tolsym*half*(norm2a+norm2b))then 2126 ! The lattice is face centered 2127 fact=2 ; center=-3 2128 cell_base(:,ia)=2.0d0*vecta(:) 2129 cell_base(:,ib)=2.0d0*vectb(:) 2130 cell_base(:,itrial)=minim(:,itrial) 2131 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2132 end if 2133 end if 2134 end if 2135 2136 ! DEBUG 2137 ! write(std_out,*)' after test_c, foundc=',foundc 2138 ! ENDDEBUG 2139 2140 ! Working hypothesis : the conventional cell is orthogonal, 2141 ! and body centered with no basis vector being an axis, 2142 ! in which case the basis vectors must be equal (even for orthorhombic) 2143 if(foundc==0 .and. orthogonal==1 .and. & 2144 & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then 2145 ! Compute the combination of the two other vectors 2146 vecta(:)=minim(:,ia)+minim(:,ib) 2147 vectb(:)=minim(:,ia)-minim(:,ib) 2148 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2149 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 2150 ! Project the trial vector on the first of the two vectors 2151 reduceda=( minim(1,itrial)*vecta(1)+ & 2152 & minim(2,itrial)*vecta(2)+ & 2153 & minim(3,itrial)*vecta(3) )/norm2a 2154 reducedb=( minim(1,itrial)*vectb(1)+ & 2155 & minim(2,itrial)*vectb(2)+ & 2156 & minim(3,itrial)*vectb(3) )/norm2b 2157 if( abs(abs(reduceda)-0.5d0)<tolsym )then 2158 ! The first vector is an axis 2159 fact=2 ; center=-1 2160 cell_base(:,ia)=vecta(:) 2161 vecta(:)=minim(:,itrial)-reduceda*vecta(:) 2162 vectb(:)=0.5d0*vectb(:) 2163 cell_base(:,ib)=vecta(:)+vectb(:) 2164 cell_base(:,itrial)=vecta(:)-vectb(:) 2165 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2166 else if( abs(abs(reducedb)-0.5d0)<tolsym )then 2167 ! The second vector is an axis 2168 fact=2 ; center=-1 2169 cell_base(:,ib)=vectb(:) 2170 vectb(:)=minim(:,itrial)-reducedb*vectb(:) 2171 vecta(:)=0.5d0*vecta(:) 2172 cell_base(:,ia)=vectb(:)+vecta(:) 2173 cell_base(:,itrial)=vectb(:)-vecta(:) 2174 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2175 end if 2176 end if 2177 2178 ! Working hypothesis : the conventional cell is orthogonal, 2179 ! and face centered, in the case where two minimal vectors are equal 2180 if(foundc==0 .and. orthogonal==1 .and. equal(itrial)==1 ) then 2181 ! Compute the combination of these two vectors 2182 vecta(:)=minim(:,ia)+minim(:,ib) 2183 vectb(:)=minim(:,ia)-minim(:,ib) 2184 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2185 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 2186 ! Project the trial vector on the two vectors 2187 reduceda=( minim(1,itrial)*vecta(1)+ & 2188 & minim(2,itrial)*vecta(2)+ & 2189 & minim(3,itrial)*vecta(3) )/norm2a 2190 reducedb=( minim(1,itrial)*vectb(1)+ & 2191 & minim(2,itrial)*vectb(2)+ & 2192 & minim(3,itrial)*vectb(3) )/norm2b 2193 if( (abs(abs(reduceda)-0.5d0)<tolsym .and. abs(reducedb)<tolsym ) .or. & 2194 & ( abs(reduceda)<tolsym .and. abs(abs(reducedb)-0.5d0)<tolsym) )then 2195 fact=2 ; center=-3 2196 cell_base(:,itrial)= & 2197 & (minim(:,itrial)-reduceda*vecta(:)-reducedb*vectb(:) )*2.0d0 2198 cell_base(:,ia)=vecta(:) 2199 cell_base(:,ib)=vectb(:) 2200 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2201 end if 2202 end if 2203 2204 ! Working hypothesis : the conventional cell is orthogonal, 2205 ! face centered, but no two vectors are on the same "square" 2206 if(foundc==0 .and. orthogonal==1)then 2207 ! Compute the combination of these two vectors 2208 vecta(:)=minim(:,ia)+minim(:,ib) 2209 vectb(:)=minim(:,ia)-minim(:,ib) 2210 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2211 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 2212 ! The trial vector length must be equal to one of these lengths 2213 if(abs(metmin(itrial,itrial)-norm2a)<tolsym*norm2a)then 2214 fact=2 ; center=-3 2215 cell_base(:,ia)=vecta(:)+minim(:,itrial) 2216 cell_base(:,ib)=vecta(:)-minim(:,itrial) 2217 ! Project vectb perpendicular to cell_base(:,ia) and cell_base(:,ib) 2218 norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2 2219 norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2 2220 reduceda=( cell_base(1,ia)*vectb(1)+ & 2221 & cell_base(2,ia)*vectb(2)+ & 2222 & cell_base(3,ia)*vectb(3) )/norm2a 2223 reducedb=( cell_base(1,ib)*vectb(1)+ & 2224 & cell_base(2,ib)*vectb(2)+ & 2225 & cell_base(3,ib)*vectb(3) )/norm2b 2226 if( abs(abs(reduceda)-0.5d0)<tolsym .and. & 2227 & abs(abs(reducedb)-0.5d0)<tolsym )then 2228 cell_base(:,itrial)=vectb(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib) 2229 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2230 end if 2231 else if(abs(metmin(itrial,itrial)-norm2b)<tolsym*norm2b)then 2232 fact=2 ; center=-3 2233 cell_base(:,ia)=vectb(:)+minim(:,itrial) 2234 cell_base(:,ib)=vectb(:)-minim(:,itrial) 2235 ! Project vecta perpendicular to cell_base(:,ia) and cell_base(:,ib) 2236 norm2a=cell_base(1,ia)**2+cell_base(2,ia)**2+cell_base(3,ia)**2 2237 norm2b=cell_base(1,ib)**2+cell_base(2,ib)**2+cell_base(3,ib)**2 2238 reduceda=( cell_base(1,ia)*vecta(1)+ & 2239 & cell_base(2,ia)*vecta(2)+ & 2240 & cell_base(3,ia)*vecta(3) )/norm2a 2241 reducedb=( cell_base(1,ib)*vecta(1)+ & 2242 & cell_base(2,ib)*vecta(2)+ & 2243 & cell_base(3,ib)*vecta(3) )/norm2b 2244 if( abs(abs(reduceda)-0.5d0)<tolsym .and. & 2245 & abs(abs(reducedb)-0.5d0)<tolsym )then 2246 cell_base(:,itrial)=vecta(:)-reduceda*cell_base(:,ia)-reducedb*cell_base(:,ib) 2247 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2248 end if 2249 end if 2250 end if 2251 2252 ! Working hypothesis : the cell is rhombohedral, and 2253 ! the three minimal vectors have same length and same absolute scalar product 2254 if(foundc==0 .and. iholohedry==5 .and. & 2255 & equal(1)==1 .and. equal(2)==1 .and. equal(3)==1 )then 2256 if( abs(abs(metmin(1,2))-abs(metmin(1,3)))<tolsym*metmin(1,1) .and. & 2257 & abs(abs(metmin(1,2))-abs(metmin(2,3)))<tolsym*metmin(2,2) )then 2258 fact=1 ; center=0 2259 cell_base(:,:)=minim(:,:) 2260 ! One might have to change the sign of one of the vectors 2261 sign12=1 ; sign13=1 ; sign23=1 2262 if(metmin(1,2)<0.0d0)sign12=-1 2263 if(metmin(1,3)<0.0d0)sign13=-1 2264 if(metmin(2,3)<0.0d0)sign23=-1 2265 sumsign=sign12+sign13+sign23 2266 if(sumsign==-1)then 2267 if(sign12==1)cell_base(:,3)=-cell_base(:,3) 2268 if(sign13==1)cell_base(:,2)=-cell_base(:,2) 2269 if(sign23==1)cell_base(:,1)=-cell_base(:,1) 2270 else if(sumsign==1)then 2271 if(sign12==-1)cell_base(:,3)=-cell_base(:,3) 2272 if(sign13==-1)cell_base(:,2)=-cell_base(:,2) 2273 if(sign23==-1)cell_base(:,1)=-cell_base(:,1) 2274 end if 2275 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2276 end if 2277 end if 2278 2279 ! DEBUG 2280 ! write(std_out,*)' after test_3a, foundc=',foundc 2281 ! write(std_out,*)' after test_3a, itrial=',itrial 2282 ! write(std_out,*)' after test_3a, equal(:)=',equal(:) 2283 ! ENDDEBUG 2284 2285 ! Working hypothesis : the cell is rhombohedral, one vector 2286 ! is parallel to the trigonal axis 2287 if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 )then 2288 vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib) 2289 norm2trial=minim(1,itrial)**2+minim(2,itrial)**2+minim(3,itrial)**2 2290 reduceda=( minim(1,itrial)*vecta(1)+ & 2291 & minim(2,itrial)*vecta(2)+ & 2292 & minim(3,itrial)*vecta(3) )/norm2trial 2293 reducedb=( minim(1,itrial)*vectb(1)+ & 2294 & minim(2,itrial)*vectb(2)+ & 2295 & minim(3,itrial)*vectb(3) )/norm2trial 2296 ! DEBUG 2297 ! write(std_out,*)' reduceda,reducedb=',reduceda,reducedb 2298 ! ENDDEBUG 2299 if(abs(abs(reduceda)-1.0d0/3.0d0)<tolsym .and. & 2300 & abs(abs(reducedb)-1.0d0/3.0d0)<tolsym ) then 2301 ! Possibly change of sign to make positive the scalar product with 2302 ! the vector parallel to the trigonal axis 2303 if(reduceda<zero)vecta(:)=-vecta(:) 2304 if(reducedb<zero)vectb(:)=-vectb(:) 2305 ! Projection on the orthogonal plane 2306 vecta(:)=vecta(:)-abs(reduceda)*cell_base(:,itrial) 2307 vectb(:)=vectb(:)-abs(reducedb)*cell_base(:,itrial) 2308 ! These two vectors should have an angle of 120 degrees 2309 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2310 scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3) 2311 ! DEBUG 2312 ! write(std_out,*)' norm2a,scalarprod=',norm2a,scalarprod 2313 ! ENDDEBUG 2314 if(abs(two*scalarprod+norm2a)<tolsym*norm2a)then 2315 fact=1 ; center=0 2316 if(scalarprod>0.0d0)vectb(:)=-vectb(:) 2317 ! Now vecta and vectb have an angle of 120 degrees 2318 cell_base(:,1)=cell_base(:,itrial)/3.0d0+vecta(:) 2319 cell_base(:,2)=cell_base(:,itrial)/3.0d0+vectb(:) 2320 cell_base(:,3)=cell_base(:,itrial)/3.0d0-vecta(:)-vectb(:) 2321 ! DEBUG 2322 ! write(std_out,*)' cell_base(:,1)=',cell_base(:,1) 2323 ! write(std_out,*)' cell_base(:,2)=',cell_base(:,2) 2324 ! write(std_out,*)' cell_base(:,3)=',cell_base(:,3) 2325 ! ENDDEBUG 2326 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2327 end if 2328 end if 2329 end if 2330 2331 ! Working hypothesis : the cell is rhombohedral, one vector 2332 ! is in the plane perpendicular to the trigonal axis 2333 if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then 2334 vecta(:)=minim(:,ia)+minim(:,ib) 2335 vectb(:)=minim(:,ia)-minim(:,ib) 2336 norm2trial=cell_base(1,itrial)**2+cell_base(2,itrial)**2+cell_base(3,itrial)**2 2337 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2338 norm2b=vecta(1)**2+vecta(2)**2+vecta(3)**2 2339 reduceda=( cell_base(1,itrial)*vecta(1)+ & 2340 & cell_base(2,itrial)*vecta(2)+ & 2341 & cell_base(3,itrial)*vecta(3) )/norm2trial 2342 reducedb=( cell_base(1,itrial)*vectb(1)+ & 2343 & cell_base(2,itrial)*vectb(2)+ & 2344 & cell_base(3,itrial)*vectb(3) )/norm2trial 2345 if(abs(norm2trial-norm2a)<tolsym*norm2a .and. & 2346 & abs(abs(2*reduceda)-norm2trial)<tolsym*norm2trial )then 2347 fact=1 ; center=0 2348 cell_base(:,1)=minim(:,ia) 2349 cell_base(:,2)=-minim(:,ib) 2350 cell_base(:,3)=-minim(:,ib)+2*reduceda*minim(:,itrial) 2351 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2352 else if (abs(norm2trial-norm2b)<tolsym*norm2b .and. & 2353 & abs(abs(2*reducedb)-norm2trial)<tolsym*norm2trial )then 2354 fact=1 ; center=0 2355 cell_base(:,1)=minim(:,ia) 2356 cell_base(:,2)=minim(:,ib) 2357 cell_base(:,3)=minim(:,ib)+2*reducedb*minim(:,itrial) 2358 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2359 end if 2360 end if 2361 2362 ! Working hypothesis : the cell is rhombohedral, two vectors 2363 ! are in the plane perpendicular to the trigonal axis 2364 if(foundc==0 .and. iholohedry==5 .and. equal(itrial)==1 ) then 2365 vecta(:)=minim(:,ia) ; vectb(:)=minim(:,ib) 2366 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2367 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 2368 scalarprod=vecta(1)*vectb(1)+vecta(2)*vectb(2)+vecta(3)*vectb(3) 2369 if(abs(abs(2*scalarprod)-norm2a)<tolsym*norm2a)then 2370 ! This is in order to have 120 angle between vecta and vectb 2371 if(scalarprod>0.0d0)vectb(:)=-vectb(:) 2372 reduceda=( cell_base(1,itrial)*vecta(1)+ & 2373 & cell_base(2,itrial)*vecta(2)+ & 2374 & cell_base(3,itrial)*vecta(3) )/norm2a 2375 reducedb=( cell_base(1,itrial)*vectb(1)+ & 2376 & cell_base(2,itrial)*vectb(2)+ & 2377 & cell_base(3,itrial)*vectb(3) )/norm2b 2378 fact=1 ; center=0 2379 cell_base(:,1)=minim(:,itrial) 2380 if(abs(reduceda-0.5d0)<tolsym .and. abs(reducedb)<tolsym )then 2381 cell_base(:,2)=minim(:,itrial)-vecta(:) 2382 cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:) 2383 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2384 else if(abs(reduceda-0.5d0)<tolsym.and. abs(reducedb+0.5d0)<tolsym )then 2385 cell_base(:,2)=minim(:,itrial)-vecta(:) 2386 cell_base(:,3)=minim(:,itrial)+vectb(:) 2387 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2388 else if(abs(reduceda)<tolsym .and. abs(reducedb+0.5d0)<tolsym )then 2389 cell_base(:,2)=minim(:,itrial)+vectb(:) 2390 cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:) 2391 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2392 else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb)<tolsym)then 2393 cell_base(:,2)=minim(:,itrial)+vecta(:) 2394 cell_base(:,3)=minim(:,itrial)+vecta(:)+vectb(:) 2395 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2396 else if(abs(reduceda+0.5d0)<tolsym .and. abs(reducedb-0.5d0)<tolsym)then 2397 cell_base(:,2)=minim(:,itrial)+vecta(:) 2398 cell_base(:,3)=minim(:,itrial)-vectb(:) 2399 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2400 else if(abs(reduceda)<tolsym .and. abs(reducedb-0.5d0)<tolsym )then 2401 cell_base(:,2)=minim(:,itrial)-vectb(:) 2402 cell_base(:,3)=minim(:,itrial)-vecta(:)-vectb(:) 2403 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2404 end if 2405 end if 2406 end if 2407 2408 ! Working hypothesis : monoclinic holohedry, primitive. Then, two angles are 90 degrees 2409 if(foundc==0 .and. iholohedry==2 .and. & 2410 & ang90(ia)==1 .and. ang90(ib)==1 ) then 2411 fact=1 ; center=0 2412 cell_base(:,1)=minim(:,ia) 2413 cell_base(:,2)=minim(:,itrial) 2414 cell_base(:,3)=minim(:,ib) 2415 ! Checks that the basis vectors are OK for the target holohedry 2416 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2417 end if 2418 2419 ! Monoclinic holohedry, one-face-centered cell 2420 ! Working hypothesis, two vectors have equal length. 2421 do icase=1,5 2422 if(foundc==0 .and. iholohedry==2 .and. equal(itrial)==1 ) then 2423 vecta(:)=cell_base(:,ia)+cell_base(:,ib) 2424 vectb(:)=cell_base(:,ia)-cell_base(:,ib) 2425 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2426 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 2427 ! The minim(:,trial) vector belongs to the 2428 ! plane parallel to the cell_base(:,ia),cell_base(:,ib) plane 2429 ! In that plane, must try minim(:,itrial), 2430 ! as well as the 4 different combinations of 2431 ! minim(:,itrial) with the vectors in the plane 2432 if(icase==1)vectc(:)=minim(:,itrial) 2433 if(icase==2)vectc(:)=minim(:,itrial)+cell_base(:,ia) 2434 if(icase==3)vectc(:)=minim(:,itrial)+cell_base(:,ib) 2435 if(icase==4)vectc(:)=minim(:,itrial)-cell_base(:,ia) 2436 if(icase==5)vectc(:)=minim(:,itrial)-cell_base(:,ib) 2437 norm2c=vectc(1)**2+vectc(2)**2+vectc(3)**2 2438 sca=vectc(1)*vecta(1)+& 2439 & vectc(2)*vecta(2)+& 2440 & vectc(3)*vecta(3) 2441 scb=vectc(1)*vectb(1)+& 2442 & vectc(2)*vectb(2)+& 2443 & vectc(3)*vectb(3) 2444 ! DEBUG 2445 ! write(std_out,*)' symlatt : test iholohedry=2, sca,scb=',sca,scb 2446 ! ENDDEBUG 2447 if(abs(sca)<tolsym*sqrt(norm2c*norm2a) .or. abs(scb)<tolsym*sqrt(norm2c*norm2b))then 2448 fact=2 ; center=3 2449 ! The itrial direction is centered 2450 cell_base(:,3)=vectc(:) 2451 if(abs(sca)<tolsym*sqrt(norm2c*norm2a))then 2452 cell_base(:,2)=vecta(:) 2453 cell_base(:,1)=vectb(:) 2454 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2455 else if(abs(scb)<tolsym*sqrt(norm2c*norm2b))then 2456 cell_base(:,2)=vectb(:) 2457 cell_base(:,1)=vecta(:) 2458 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2459 end if 2460 end if 2461 end if 2462 end do ! icase=1,5 2463 2464 ! Monoclinic holohedry, one-face-centered cell, but non equivalent. 2465 ! This case, one pair of vectors is orthogonal 2466 if(foundc==0 .and. iholohedry==2 .and. ang90(itrial)==1) then 2467 vecta(:)=minim(:,ia) 2468 vectb(:)=minim(:,ib) 2469 norm2a=vecta(1)**2+vecta(2)**2+vecta(3)**2 2470 norm2b=vectb(1)**2+vectb(2)**2+vectb(3)**2 2471 ! Project the trial vector on the two vectors 2472 reduceda=( minim(1,itrial)*vecta(1)+ & 2473 & minim(2,itrial)*vecta(2)+ & 2474 & minim(3,itrial)*vecta(3) )/norm2a 2475 reducedb=( minim(1,itrial)*vectb(1)+ & 2476 & minim(2,itrial)*vectb(2)+ & 2477 & minim(3,itrial)*vectb(3) )/norm2b 2478 if(abs(abs(reduceda)-0.5d0)<tolsym .or. abs(abs(reducedb)-0.5d0)<tolsym) then 2479 fact=2 ; center=3 2480 if(abs(abs(reduceda)-0.5d0)<tolsym)then 2481 cell_base(:,2)=vecta(:) 2482 cell_base(:,3)=vectb(:) 2483 cell_base(:,1)=2*(minim(:,itrial)-reduceda*vecta(:)) 2484 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2485 else if(abs(abs(reducedb)-0.5d0)<tolsym)then 2486 cell_base(:,2)=vectb(:) 2487 cell_base(:,3)=vecta(:) 2488 cell_base(:,1)=2*(minim(:,itrial)-reducedb*vectb(:)) 2489 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2490 end if 2491 end if 2492 end if 2493 2494 ! Monoclinic holohedry, one-face-centered cell, but non equivalent. 2495 ! This case, no pair of vectors is orthogonal, no pair of vector of equal lentgh 2496 if(foundc==0 .and. iholohedry==2)then 2497 ! Try to find a vector that belongs to the mediator plane, or the binary vector. 2498 ! There must be one such vector, if centered monoclinic and no pair of vectors of equal length, 2499 ! either among the three vectors, or among one of their differences or sums. 2500 ! And there must be, among the two other vectors, one vector whose projection 2501 ! on this vector is half the length of this vector. 2502 vecta(:)=minim(:,ia) 2503 vectb(:)=minim(:,ib) 2504 ! Try the different possibilities for the vector on which the projection will be half ... 2505 do ii=1,5 2506 if(ii==1)vectc(:)=minim(:,itrial) 2507 if(ii==2)vectc(:)=minim(:,itrial)+vecta(:) 2508 if(ii==3)vectc(:)=minim(:,itrial)-vecta(:) 2509 if(ii==4)vectc(:)=minim(:,itrial)+vectb(:) 2510 if(ii==5)vectc(:)=minim(:,itrial)-vectb(:) 2511 norm2trial=vectc(1)**2+vectc(2)**2+vectc(3)**2 2512 ! Project the two vectors on the trial vector 2513 reduceda=( vectc(1)*vecta(1)+vectc(2)*vecta(2)+vectc(3)*vecta(3) )/norm2trial 2514 reducedb=( vectc(1)*vectb(1)+vectc(2)*vectb(2)+vectc(3)*vectb(3) )/norm2trial 2515 found=0 2516 if(abs(abs(reduceda)-0.5d0)<tolsym)then 2517 vin1(:)=vectc(:) 2518 vin2(:)=2.0d0*(vecta(:)-reduceda*vectc(:)) 2519 vext(:)=vectb(:) 2520 found=1 2521 else if(abs(abs(reducedb)-0.5d0)<tolsym)then 2522 vin1(:)=vectc(:) 2523 vin2(:)=2.0d0*(vectb(:)-reduceda*vectc(:)) 2524 vext(:)=vecta(:) 2525 found=1 2526 end if 2527 if(found==1)exit 2528 end do 2529 ! Now, vin1 and vin2 are perpendicular to each other, and in the plane that contains the binary vector. 2530 ! One of them must be the binary vector if any. 2531 ! On the other hand, vext is out-of-plane. Might belong to the mediator plane or not. 2532 ! If C monoclinc, then the projection of this vext on the binary vector will be either 0 or +1/2 or -1/2. 2533 ! The binary axis must be stored in cell_base(:,2) for conventional C-cell 2534 if(found==1)then 2535 found=0 2536 2537 ! Test vin1 being the binary axis 2538 norm2trial=vin1(1)**2+vin1(2)**2+vin1(3)**2 2539 reduceda=(vext(1)*vin1(1)+vext(2)*vin1(2)+vext(3)*vin1(3))/norm2trial 2540 if(abs(reduceda)<tolsym)then ! vin1 is the binary axis and vext is in the mediator plane 2541 found=1 2542 cell_base(:,1)=vin2(:) 2543 cell_base(:,2)=vin1(:) 2544 cell_base(:,3)=vext(:) 2545 else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin1 is the binary axis and vext has +1/2 or -1/2 as projection 2546 found=1 2547 cell_base(:,1)=vin2(:) 2548 cell_base(:,2)=vin1(:) 2549 cell_base(:,3)=vext(:)-reduceda*vin1(:)+vin2(:)*half 2550 else 2551 ! Test vin2 being the binary axis 2552 norm2trial=vin2(1)**2+vin2(2)**2+vin2(3)**2 2553 reduceda=(vext(1)*vin2(1)+vext(2)*vin2(2)+vext(3)*vin2(3))/norm2trial 2554 if(abs(reduceda)<tolsym)then ! vin2 is the binary axis and vext is in the mediator plane 2555 found=1 2556 cell_base(:,1)=vin1(:) 2557 cell_base(:,2)=vin2(:) 2558 cell_base(:,3)=vext(:) 2559 else if(abs(abs(reduceda)-0.5d0)<tolsym)then ! vin2 is the binary axis and vext has +1/2 or -1/2 as projection 2560 found=1 2561 cell_base(:,1)=vin1(:) 2562 cell_base(:,2)=vin2(:) 2563 cell_base(:,3)=vext(:)-reduceda*vin2(:)+vin1(:)*half 2564 end if 2565 end if 2566 2567 if(found==1)then 2568 fact=2 ; center=3 2569 call holocell(cell_base,0,foundc,iholohedry,tolsym) 2570 end if 2571 end if 2572 end if 2573 2574 end do ! Do-loop on three different directions 2575 end do ! Do-loop on different target holohedries 2576 2577 if(foundc==0)then 2578 iholohedry=1 ; fact=1 ; center=0 2579 cell_base(:,:)=minim(:,:) 2580 end if 2581 2582 !DEBUG 2583 !write(std_out,*)' symlatt : done with centering tests, foundc=',foundc 2584 !write(std_out,*)' center=',center 2585 !write(std_out,*)' iholohedry=',iholohedry 2586 !call flush(std_out) 2587 !ENDDEBUG 2588 2589 !-------------------------------------------------------------------------- 2590 !Final check on the Bravais lattice, using the basis vectors 2591 2592 !Recompute the metric tensor 2593 if(foundc==1)then 2594 do ii=1,3 2595 metmin(:,ii)=cell_base(1,:)*cell_base(1,ii)+& 2596 & cell_base(2,:)*cell_base(2,ii)+& 2597 & cell_base(3,:)*cell_base(3,ii) 2598 end do 2599 end if 2600 2601 !Examine the angles and vector lengths 2602 ang90(:)=0 2603 if(metmin(1,2)**2<tolsym**2*metmin(1,1)*metmin(2,2))ang90(3)=1 2604 if(metmin(1,3)**2<tolsym**2*metmin(1,1)*metmin(3,3))ang90(2)=1 2605 if(metmin(2,3)**2<tolsym**2*metmin(2,2)*metmin(3,3))ang90(1)=1 2606 equal(:)=0 2607 if(abs(metmin(1,1)-metmin(2,2))<tolsym*half*(metmin(1,1)+metmin(2,2)))equal(3)=1 2608 if(abs(metmin(1,1)-metmin(3,3))<tolsym*half*(metmin(1,1)+metmin(3,3)))equal(2)=1 2609 if(abs(metmin(2,2)-metmin(3,3))<tolsym*half*(metmin(2,2)+metmin(3,3)))equal(1)=1 2610 2611 !DEBUG 2612 !write(std_out,*)' symlatt : recompute the metric tensor ' 2613 !write(std_out,*)' ang90=',ang90 2614 !write(std_out,*)' equal=',equal 2615 !call flush(std_out) 2616 !ENDDEBUG 2617 2618 !The axes will be aligned with the previously determined 2619 !basis vectors, EXCEPT for the tetragonal cell, see later 2620 axes(:,:)=cell_base(:,:) 2621 2622 found=0 2623 !Check orthogonal conventional cells 2624 if(ang90(1)+ang90(2)+ang90(3)==3)then 2625 2626 ! Cubic system 2627 if(equal(1)+equal(2)+equal(3)==3)then 2628 ! However, one-face centered is not admitted 2629 if(center==0 .or. center==-1 .or. center==-3)then 2630 iholohedry=7 ; found=1 2631 if(center==0)then 2632 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cP (primitive cubic)' 2633 else if(center==-1)then 2634 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cI (body-centered cubic)' 2635 else if(center==-3)then 2636 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is cF (face-centered cubic)' 2637 end if 2638 end if 2639 end if 2640 2641 ! Tetragonal system 2642 if(found==0 .and. & 2643 & (equal(1)==1 .or. equal(2)==1 .or. equal(3)==1) )then 2644 ! However, one-face centered or face-centered is not admitted 2645 if(center==0 .or. center==-1)then 2646 iholohedry=4 ; found=1 2647 if(equal(1)==1)then 2648 axes(:,3)=cell_base(:,1) ; axes(:,1)=cell_base(:,2) ; axes(:,2)=cell_base(:,3) 2649 else if(equal(2)==1)then 2650 axes(:,3)=cell_base(:,2) ; axes(:,2)=cell_base(:,1) ; axes(:,1)=cell_base(:,3) 2651 else if(equal(3)==1)then 2652 axes(:,:)=cell_base(:,:) 2653 end if 2654 if(center==0)then 2655 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tP (primitive tetragonal)' 2656 else if(center==-1)then 2657 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is tI (body-centered tetragonal)' 2658 end if 2659 end if 2660 end if 2661 2662 ! Orthorhombic system 2663 if(found==0)then 2664 iholohedry=3 ; found=1 2665 axes(:,:)=cell_base(:,:) 2666 if(center==0)then 2667 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oP (primitive orthorhombic)' 2668 else if(center==-1)then 2669 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oI (body-centered orthorhombic)' 2670 else if(center==1 .or. center==2 .or. center==3)then 2671 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oC (one-face-centered orthorhombic)' 2672 else if(center==-3)then 2673 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is oF (face-centered orthorhombic)' 2674 end if 2675 end if 2676 2677 else 2678 2679 ! Hexagonal system 2680 if(found==0 .and. ang90(1)==1 .and. ang90(2)==1 .and. equal(3)==1 .and. (2*metmin(2,1)+metmin(1,1))<tolsym*metmin(1,1))then 2681 iholohedry=6 ; found=1 2682 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hP (primitive hexagonal)' 2683 end if 2684 2685 ! Rhombohedral system 2686 if(found==0 .and. equal(1)+equal(2)+equal(3)==3 .and. & 2687 & abs(metmin(2,1)-metmin(3,2))<tolsym*metmin(2,2) .and. & 2688 & abs(metmin(2,1)-metmin(3,1))<tolsym*metmin(1,1) )then 2689 iholohedry=5 ; found=1 2690 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is hR (rhombohedral)' 2691 end if 2692 2693 ! Monoclinic system 2694 if(found==0 .and. ang90(1)+ang90(2)+ang90(3)==2 )then 2695 iholohedry=2 ; found=1 2696 if(center==0)then 2697 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mP (primitive monoclinic)' 2698 else if(center==3)then 2699 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is mC (one-face-centered monoclinic)' 2700 end if 2701 end if 2702 2703 ! Triclinic system 2704 if(found==0)then 2705 iholohedry=1 ; found=1 2706 write(msg,'(a,a)')ch10,' symlatt: the Bravais lattice is aP (primitive triclinic)' 2707 end if 2708 2709 end if 2710 2711 call wrtout(iout,msg) 2712 2713 !DEBUG 2714 !write(std_out,*)' symlatt : after checking conventional orthogonal cell ' 2715 !call flush(std_out) 2716 !ENDDEBUG 2717 2718 !-------------------------------------------------------------------------- 2719 !Make sure that axes form a right-handed coordinate system 2720 !(Note : this should be done in the body of the routine, 2721 !by making changes that leave the sign of the mixed product of the three 2722 !vectors invariant) 2723 determinant=axes(1,1)*axes(2,2)*axes(3,3) & 2724 & +axes(1,2)*axes(2,3)*axes(3,1) & 2725 & +axes(1,3)*axes(3,2)*axes(2,1) & 2726 & -axes(1,1)*axes(3,2)*axes(2,3) & 2727 & -axes(1,3)*axes(2,2)*axes(3,1) & 2728 & -axes(1,2)*axes(2,1)*axes(3,3) 2729 if(determinant<zero)then 2730 axes(:,:)=-axes(:,:) 2731 end if 2732 2733 !DEBUG 2734 !write(std_out,'(a,i4)')' symlatt : before itrial do loop, iholohedry= ',iholohedry 2735 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',& 2736 !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3) 2737 !call flush(std_out) 2738 !ENDDEBUG 2739 2740 !-------------------------------------------------------------------------- 2741 !Prefer symmetry axes on the same side as the primitive axes, 2742 !when the changes are allowed 2743 do itrial=1,100 2744 2745 ! DEBUG 2746 ! write(std_out,'(a)')' ' 2747 ! write(std_out,'(a,i5)')' symlatt : itrial do loop, itrial= ',itrial 2748 ! write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',& 2749 ! & axes(:,1),ch10,axes(:,2),ch10,axes(:,3) 2750 ! call flush(std_out) 2751 ! ENDDEBUG 2752 2753 do ia=1,3 2754 scprods(ia,:)=axes(1,ia)*rprimd(1,:)+& 2755 & axes(2,ia)*rprimd(2,:)+& 2756 & axes(3,ia)*rprimd(3,:) 2757 norm2trial=sum(axes(:,ia)**2) 2758 scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial) 2759 end do 2760 do ia=1,3 2761 norm2trial=sum(rprimd(:,ia)**2) 2762 scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial) 2763 end do 2764 2765 !DEBUG 2766 ! write(std_out,'(a,3f12.6)')' diagonal scalar products ',scprods(1,1),scprods(2,2),scprods(3,3) 2767 ! call flush(std_out) 2768 !ENDDEBUG 2769 2770 ! One should now try all the generators of the 2771 ! proper rotations of each Bravais lattice, coupled with change of 2772 ! signs of each vector. This is not done systematically in what follows ... 2773 ! Here, the third axis is left unchanged 2774 if(iholohedry/=5)then 2775 if(scprods(1,1)<-tolsym .and. scprods(2,2)<-tolsym)then 2776 axes(:,1)=-axes(:,1) ; axes(:,2)=-axes(:,2) 2777 cycle 2778 end if 2779 end if 2780 ! The first (or second) axis is left unchanged 2781 if(iholohedry/=5 .and. iholohedry/=6)then 2782 if(scprods(2,2)<-tolsym .and. scprods(3,3)<-tolsym)then 2783 axes(:,2)=-axes(:,2) ; axes(:,3)=-axes(:,3) 2784 cycle 2785 end if 2786 if(scprods(1,1)<-tolsym .and. scprods(3,3)<-tolsym)then 2787 axes(:,1)=-axes(:,1) ; axes(:,3)=-axes(:,3) 2788 cycle 2789 end if 2790 end if 2791 ! Permutation of the three axis 2792 if(iholohedry==5 .or. iholohedry==7)then 2793 trace=scprods(1,1)+scprods(2,2)+scprods(3,3) 2794 if(trace+tolsym< scprods(1,2)+scprods(2,3)+scprods(3,1))then 2795 vecta(:)=axes(:,1) ; axes(:,1)=axes(:,3) 2796 axes(:,3)=axes(:,2); axes(:,2)=vecta(:) 2797 cycle 2798 end if 2799 if(trace+tolsym < scprods(1,3)+scprods(2,1)+scprods(3,2))then 2800 vecta(:)=axes(:,1) ; axes(:,1)=axes(:,2) 2801 axes(:,2)=axes(:,3); axes(:,3)=vecta(:) 2802 cycle 2803 end if 2804 ! This case is observed when the three new vectors 2805 ! are pointing opposite to the three original vectors 2806 ! One takes their opposite, then switch two of them, then process 2807 ! them again in the loop 2808 if(sum(scprods(:,:))<-tolsym)then 2809 axes(:,1)=-axes(:,1) 2810 vecta(:)=-axes(:,2) 2811 axes(:,2)=-axes(:,3) 2812 axes(:,3)=vecta(:) 2813 cycle 2814 end if 2815 end if 2816 2817 ! Actually, for iholohedry==7, can test specifically all possibilities 2818 ! and take the best one. 2819 ! Not activated, because changing the order of symmetries in many tests ! 2820 if(iholohedry==7 .and. .false.)then 2821 ! if(iholohedry==7)then 2822 2823 !DEBUG 2824 !write(std_out,'(a,a)')ch10,' enter search of all possibilities for iholohedry==7 ' 2825 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',& 2826 !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3) 2827 !call flush(std_out) 2828 !ENDDEBUG 2829 2830 do iaxis1=1,3 2831 axes_try(:,1)=axes(:,iaxis1) 2832 do iaxis2=1,2 2833 if(iaxis1==1)axes_try(:,2)=axes(:,1+iaxis2) 2834 if(iaxis1==2)axes_try(:,2)=axes(:,2*iaxis2-1) 2835 if(iaxis1==3)axes_try(:,2)=axes(:,iaxis2) 2836 if(iaxis2==1)axes_try(:,3)=axes(:,3) 2837 if(iaxis2==2)axes_try(:,3)=axes(:,1) 2838 if(iaxis1==1.and.iaxis2==2)axes_try(:,3)=axes(:,2) 2839 if(iaxis1==3.and.iaxis2==1)axes_try(:,3)=axes(:,2) 2840 do isign1=1,-1,-2 2841 axes_try(:,1)=-axes_try(:,1) 2842 do isign2=1,-1,-2 2843 axes_try(:,2)=-axes_try(:,2) 2844 determinant=axes_try(1,1)*axes_try(2,2)*axes_try(3,3) & 2845 & +axes_try(1,2)*axes_try(2,3)*axes_try(3,1) & 2846 & +axes_try(1,3)*axes_try(3,2)*axes_try(2,1) & 2847 & -axes_try(1,1)*axes_try(3,2)*axes_try(2,3) & 2848 & -axes_try(1,3)*axes_try(2,2)*axes_try(3,1) & 2849 & -axes_try(1,2)*axes_try(2,1)*axes_try(3,3) 2850 if(determinant<zero)axes_try(:,3)=-axes_try(:,3) 2851 do ia=1,3 2852 scprods(ia,:)=axes_try(1,ia)*rprimd(1,:)+& 2853 & axes_try(2,ia)*rprimd(2,:)+& 2854 & axes_try(3,ia)*rprimd(3,:) 2855 norm2trial=sum(axes_try(:,ia)**2) 2856 scprods(ia,:)=scprods(ia,:)/sqrt(norm2trial) 2857 end do 2858 do ia=1,3 2859 norm2trial=sum(rprimd(:,ia)**2) 2860 scprods(:,ia)=scprods(:,ia)/sqrt(norm2trial) 2861 end do 2862 trace=scprods(1,1)+scprods(2,2)+scprods(3,3) 2863 if(iaxis1==1.and.iaxis2==1.and.isign1==1.and.isign2==1)then 2864 trace_best=trace 2865 axes_best=axes_try 2866 else if (trace>trace_best+tolsym)then 2867 trace_best=trace 2868 axes_best=axes_try 2869 endif 2870 enddo ! isign2 2871 enddo ! isign1 2872 enddo ! iaxes2 2873 enddo ! iaxes1 2874 axes=axes_best 2875 endif ! iholohedry=7 2876 exit 2877 end do 2878 2879 !-------------------------------------------------------------------------- 2880 2881 !DEBUG 2882 !write(std_out,'(a,a)')ch10,' after order/sign optimization do-loop ' 2883 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',& 2884 !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3) 2885 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' axes =',& 2886 !& axes(:,1),ch10,axes(:,2),ch10,axes(:,3) 2887 !call flush(std_out) 2888 !ENDDEBUG 2889 2890 !Compute the coordinates of rprimd in the system defined by axes(:,:) 2891 call matr3inv(axes,axesinvt) 2892 do ii=1,3 2893 coord(:,ii)=rprimd(1,ii)*axesinvt(1,:)+ & 2894 & rprimd(2,ii)*axesinvt(2,:)+ & 2895 & rprimd(3,ii)*axesinvt(3,:) 2896 end do 2897 2898 !Check that the coordinates are integers, or half-integer in 2899 !the case there is a centering, and generate integer coordinates 2900 do ii=1,3 2901 do jj=1,3 2902 val=coord(ii,jj)*fact 2903 if(abs(val-nint(val))>fact*two*tolsym)then 2904 write(msg,'(4a,a,3es18.10,a,a,3es18.10,a,a,3es18.10,a,a,i4)')& 2905 'One of the coordinates of rprimd in axes is non-integer,',ch10,& 2906 'or non-half-integer (if centering), within 2*tolsym.',ch10,& 2907 'coord=',coord(:,1),ch10,& 2908 ' ',coord(:,2),ch10,& 2909 ' ',coord(:,3),ch10,& 2910 'fact=',fact 2911 ABI_BUG(msg) 2912 end if 2913 icoord(ii,jj)=nint(val) 2914 end do 2915 end do 2916 2917 !Store the bravais lattice characteristics 2918 bravais(1)=iholohedry 2919 bravais(2)=center 2920 bravais(3:5)=icoord(1:3,1) 2921 bravais(6:8)=icoord(1:3,2) 2922 bravais(9:11)=icoord(1:3,3) 2923 2924 !-------------------------------------------------------------------------- 2925 !Initialize the set of symmetries 2926 !Bravais lattices are always invariant under identity and inversion 2927 2928 !Identity and inversion 2929 ptsymrel(:,:,1)=identity(:,:) ; ptsymrel(:,:,2)=-identity(:,:) 2930 nptsym=2 2931 2932 !Keep this for IFCv70 compiler 2933 if(nptsym/=2)then 2934 write(msg,'(a,a,a,a)')ch10,& 2935 ' symlatt : BUG -',ch10,& 2936 ' Crazy error, compiler bug ' 2937 call wrtout(std_out,msg) 2938 end if 2939 2940 !-------------------------------------------------------------------------- 2941 !Initialize some generators 2942 !gen6 is defined in a coordinated system with gamma=120 degrees 2943 gen6(:,:)=0 ; gen6(3,3)=1 ; gen6(1,1)=1 ; gen6(1,2)=-1 ; gen6(2,1)=1 2944 gen3(:,:)=0 ; gen3(1,2)=1 ; gen3(2,3)=1 ; gen3(3,1)=1 2945 gen2xy(:,:)=0 ; gen2xy(2,1)=1 ; gen2xy(1,2)=1; gen2xy(3,3)=1 2946 gen2y(:,:)=0 ; gen2y(1,1)=-1; gen2y(2,2)=1 ; gen2y(3,3)=-1 2947 gen2z(:,:)=0 ; gen2z(1,1)=-1; gen2z(2,2)=-1; gen2z(3,3)=1 2948 2949 !-------------------------------------------------------------------------- 2950 2951 !Define the generators for each holohedry (inversion is already included) 2952 if(iholohedry==6)then 2953 ngen=2 2954 gen(:,:,1)=gen2xy(:,:) ; order(1)=2 2955 gen(:,:,2)=gen6(:,:) ; order(2)=6 2956 else if(iholohedry==5)then 2957 ngen=2 2958 gen(:,:,1)=gen2xy(:,:) ; order(1)=2 2959 gen(:,:,2)=gen3(:,:) ; order(2)=3 2960 else 2961 gen(:,:,1)=gen2y(:,:) ; order(1)=2 2962 gen(:,:,2)=gen2z(:,:) ; order(2)=2 2963 gen(:,:,3)=gen2xy(:,:) ; order(3)=2 2964 gen(:,:,4)=gen3(:,:) ; order(4)=3 2965 if(iholohedry<=4)ngen=iholohedry-1 2966 if(iholohedry==7)ngen=4 2967 end if 2968 2969 !Build the point symmetry operations from generators, in the reduced system 2970 !of coordinates defined by axes(:,:) 2971 if(ngen/=0)then 2972 do igen=1,ngen 2973 do isym=1+nptsym,order(igen)*nptsym 2974 jsym=isym-nptsym 2975 do ii=1,3 2976 ptsymrel(:,ii,isym)=gen(:,1,igen)*ptsymrel(1,ii,jsym)+ & 2977 & gen(:,2,igen)*ptsymrel(2,ii,jsym)+ & 2978 & gen(:,3,igen)*ptsymrel(3,ii,jsym) 2979 end do 2980 end do 2981 nptsym=order(igen)*nptsym 2982 2983 end do 2984 end if 2985 2986 !-------------------------------------------------------------------------- 2987 2988 !Transform symmetry matrices in the system defined by rprimd 2989 call symrelrot(nptsym,axes,rprimd,ptsymrel,tolsym) 2990 2991 !DEBUG 2992 !write(std_out,'(a)') ' symlatt : exit ' 2993 !call flush(std_out) 2994 !stop 2995 !ENDDEBUG 2996 2997 end subroutine symlatt
m_symfind/symspgr [ Functions ]
[ Top ] [ m_symfind ] [ Functions ]
NAME
symspgr
FUNCTION
Find the type of each symmetry operation (calling symcharac): proper symmetries 1,2,2_1,3,3_1,3_2,4,4_1,4_2,4_3,6,6_1,...6_5 improper symmetries -1,m,a,b,c,d,n,g,-3,-4,-6 , Then, build an array with the number of such operations. Then, call symlist to identify the space group.
INPUTS
bravais(11): bravais(1)=iholohedry
bravais(2)=center
bravais(3:11)=coordinates of rprimd in the axes
of the conventional bravais lattice (*2 if center/=0)
nsym=actual number of symmetries
symrel(3,3,nsym)= nsym symmetry operations in real space in terms
of primitive translations
tnons(3,nsym)=nonsymmorphic translations for each symmetry (would
be 0 0 0 each for a symmorphic space group)
OUTPUT
labels(maxsym=192)= labels of the symmetry operations spgroup=symmetry space group number
NOTES
It is assumed that the symmetry operations will be entered in the symrel tnons arrays, for the PRIMITIVE cell. The matrix of transformation from the primitive cell to the conventional cell is described in the array "bravais" (see symlatt.F90). The present routine first make the transformation from the primitive coordinates to the conventional ones, then eventually generate additional symmetries, taking into account the centering translations. Then, the order and determinant of each symmetry operation is determined. For proper symmetries (rotations), the associated translation is also determined. However, left or right handed screw rotations are not (presently) distinguished, and will be attributed equally to left or right. For the detailed description of the labelling of the axes, see symaxes.f and symplanes.f
SOURCE
1495 subroutine symspgr(bravais,labels,nsym,spgroup,symrel,tnons,tolsym) 1496 1497 use m_numeric_tools, only : OPERATOR(.x.) 1498 1499 !Arguments ------------------------------------ 1500 !scalars 1501 integer,intent(in) :: nsym 1502 integer,intent(out) :: spgroup 1503 real(dp),intent(in) :: tolsym 1504 !arrays 1505 integer,intent(in) :: bravais(11),symrel(3,3,nsym) 1506 real(dp),intent(in) :: tnons(3,nsym) 1507 character(len=128),intent(out) :: labels(192) ! 192 = maxsym 1508 1509 !Local variables------------------------------- 1510 !scalars 1511 ! logical,parameter :: verbose=.FALSE. 1512 integer :: additional_info,brvltt,center,direction=0,found,iholohedry,ii 1513 integer :: ishift,isym,jj,nshift,nsymconv,spgaxor,spgorig,sporder 1514 character(len=1) :: brvsb 1515 character(len=15) :: intsb,ptintsb,ptschsb,schsb 1516 character(len=35) :: intsbl 1517 character(len=500) :: msg 1518 !arrays 1519 integer :: ivec1(3), ivec2(3) 1520 integer :: n_axes(31),n_axest(31),prime(5),test_direction(3),symrel_uni(3,3) 1521 integer :: uniaxis(3),uniaxis_try(3) 1522 integer,allocatable :: determinant(:),symrelconv(:,:,:),t_axes(:) 1523 real(dp) :: axes(3,3),rprimdconv(3,3),vect(3,3) 1524 real(dp),allocatable :: shift(:,:),tnonsconv(:,:) 1525 1526 !************************************************************************** 1527 1528 DBG_ENTER("COLL") 1529 1530 !Initialize brvltt, from bravais(2) and bravais(1) 1531 center=bravais(2) 1532 iholohedry=bravais(1) 1533 brvltt=1 1534 if(center==-1)brvltt=2 ! Inner centering 1535 if(center==-3)brvltt=3 ! Face centering 1536 if(center==1)brvltt=5 ! A-Face centering 1537 if(center==2)brvltt=6 ! B-Face centering 1538 if(center==3)brvltt=4 ! C-Face centering 1539 if(iholohedry==5)brvltt=7 ! Rhombohedral 1540 1541 !Produce the symmetry operations, in the axis of the conventional cell 1542 nsymconv=nsym 1543 if(center/=0)nsymconv=2*nsymconv 1544 if(center==-3)nsymconv=4*nsym 1545 ABI_MALLOC(symrelconv,(3,3,nsymconv)) 1546 ABI_MALLOC(tnonsconv,(3,nsymconv)) 1547 1548 !Produce symrel and tnons in conventional axes, 1549 !name them symrelconv and tnonsconv 1550 rprimdconv(:,1)=bravais(3:5) 1551 rprimdconv(:,2)=bravais(6:8) 1552 rprimdconv(:,3)=bravais(9:11) 1553 1554 if(center/=0)rprimdconv(:,:)=rprimdconv(:,:)*half 1555 1556 axes(:,:)=zero 1557 axes(1,1)=one ; axes(2,2)=one ; axes(3,3)=one 1558 symrelconv(:,:,1:nsym)=symrel(:,:,1:nsym) 1559 !Note that the number of symmetry operations is still nsym 1560 call symrelrot(nsym,rprimdconv,axes,symrelconv,tolsym) 1561 1562 call xred2xcart(nsym,rprimdconv,tnonsconv,tnons) 1563 !Gives the associated translation, with components in the 1564 !interval ]-0.5,0.5] . 1565 tnonsconv(:,1:nsym)=tnonsconv(:,1:nsym)-nint(tnonsconv(:,1:nsym)-tol6) 1566 1567 !If the Bravais lattice is centered, duplicate or quadruplicate 1568 !the number of symmetry operations, using the Bravais 1569 !lattice shifts 1570 nshift=1 1571 if(center/=0)nshift=2 1572 if(center==-3)nshift=4 1573 ABI_MALLOC(shift,(3,nshift)) 1574 shift(:,1)=zero 1575 if(center/=0 .and. center/=-3)then 1576 shift(:,2)=half 1577 if(center==1)shift(1,2)=zero 1578 if(center==2)shift(2,2)=zero 1579 if(center==3)shift(3,2)=zero 1580 else if(center==-3)then 1581 shift(:,2)=half ; shift(1,2)=zero 1582 shift(:,3)=half ; shift(2,3)=zero 1583 shift(:,4)=half ; shift(3,4)=zero 1584 end if ! center/=0 or -3 1585 if(nshift/=1)then 1586 do ishift=2,nshift 1587 symrelconv(:,:,(ishift-1)*nsym+1:ishift*nsym)=symrelconv(:,:,1:nsym) 1588 do isym=1,nsym 1589 tnonsconv(:,(ishift-1)*nsym+isym)=tnonsconv(:,isym)+shift(:,ishift) 1590 end do 1591 end do ! ishift 1592 end if ! nshift/=1 1593 1594 !At this stage, all the symmetry operations are available, 1595 !expressed in the conventional axis, and also include 1596 !the Bravais lattive translations, and associated operations... 1597 1598 n_axes(:)=0 1599 1600 ABI_MALLOC(determinant,(nsymconv)) 1601 1602 !Get the determinant 1603 call symdet(determinant,nsymconv,symrelconv) 1604 1605 !Get the order of each the symmetry operation, as well as the maximal order 1606 !Also, examine whether each symmetry operation is the inversion, or a root 1607 !of the inversion (like -3) 1608 !Decide which kind of point symmetry operation it is 1609 !Finally assign tnonsconv order and decide the space symmetry operation 1610 1611 ABI_MALLOC(t_axes,(nsymconv)) 1612 1613 do isym=1,nsymconv 1614 1615 ! Note : nsymconv might be bigger than 192, but only for non-primitive cells, in which case labels will not be echoed anywhere. 1616 ! 192 is the fixed dimension of labels, so this avoids possible memory problems. 1617 call symcharac(center, determinant(isym), iholohedry, isym, labels(mod(isym-1,192)+1), & 1618 symrelconv(:,:,isym), tnonsconv(:,isym), t_axes(isym)) 1619 if (t_axes(isym) == -1) then 1620 write(msg, '(a,a,i3,a,3(a,3i4,a),a,3es22.12,a,a,3es22.12)' )ch10,& 1621 ' symspgr: problem with isym=',isym,ch10,& 1622 ' symrelconv(:,1,isym)=',symrelconv(:,1,isym),ch10,& 1623 ' symrelconv(:,2,isym)=',symrelconv(:,2,isym),ch10,& 1624 ' symrelconv(:,3,isym)=',symrelconv(:,3,isym),ch10,& 1625 ' tnonsconv(:,isym)=',tnonsconv(:,isym) 1626 call wrtout(std_out,msg) 1627 write(msg, '(a,i4,2a)' )& 1628 'The space symmetry operation number',isym,ch10,'is not a (translated) root of unity' 1629 ABI_BUG(msg) 1630 else if (t_axes(isym) == -2) then 1631 write(msg, '(a,i0,a)' )'The symmetry operation number ',isym,' is not a root of unity' 1632 ABI_BUG(msg) 1633 end if 1634 1635 n_axes(t_axes(isym))=n_axes(t_axes(isym))+1 1636 1637 end do ! isym=1,nsymconv 1638 1639 if (sum(n_axes)-nsymconv/=0) then 1640 write(msg, '(7a)' )& 1641 & 'Not all the symmetries have been recognized. ',ch10,& 1642 & 'This might be due either to an error in the input file',ch10,& 1643 & 'or to a BUG in ABINIT',ch10,& 1644 & 'Please contact the ABINIT group.' 1645 ABI_WARNING(msg) 1646 end if 1647 1648 !DEBUG 1649 !write(std_out,*)' symspgr : brvltt,nsymconv=',brvltt,nsymconv 1650 !write(std_out,*)' n_axes(1:10)=',n_axes(1:10) 1651 !write(std_out,*)' n_axes(11:20)=',n_axes(11:20) 1652 !write(std_out,*)' n_axes(21:31)=',n_axes(21:31) 1653 !ENDDEBUG 1654 1655 !Treat cases in which the space group cannot be identified on the 1656 !basis of n_axes one need additional information 1657 if(brvltt==1)then 1658 ! If the bravais lattice is primitive 1659 if(nsymconv==4)then 1660 n_axest=(/0,0,0,0,0,0,0,1,1,0, 0,0,0,0,0,2,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0/) 1661 if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 27 (Pcc2) or 32 (Pba2) 1662 write(std_out,*)' symspgr: 27 or 32' 1663 additional_info=2 1664 ! Select binary axis 1665 do isym=1,nsymconv 1666 if(t_axes(isym)==8)then 1667 ! Find direction of binary axis 1668 if(symrelconv(1,1,isym)==1)direction=1 1669 if(symrelconv(2,2,isym)==1)direction=2 1670 if(symrelconv(3,3,isym)==1)direction=3 1671 end if 1672 end do 1673 ! Examine the projection of the translation vector of the a, b or c mirror planes 1674 ! onto the binary axis 1675 do isym=1,nsymconv 1676 if(t_axes(isym)==16)then 1677 if(abs(tnonsconv(direction,isym))>tol8)additional_info=1 1678 end if 1679 end do 1680 end if 1681 else if(nsymconv==8)then 1682 n_axest=(/0,0,0,0,1,0,0,1,1,0, 0,0,0,0,1,2,0,0,0,2, 0,0,0,0,0,0,0,0,0,0,0/) 1683 if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 55 (Pbam) or 57 (Pbcm) 1684 write(std_out,*)' symspgr: 55 or 57' 1685 additional_info=1 1686 ! Select mirror plane m 1687 do isym=1,nsymconv 1688 if(t_axes(isym)==15)then 1689 ! Find direction of mirror plane 1690 if(symrelconv(1,1,isym)==-1)direction=1 1691 if(symrelconv(2,2,isym)==-1)direction=2 1692 if(symrelconv(3,3,isym)==-1)direction=3 1693 end if 1694 end do 1695 ! Examine the projection of the translation vector of the a, b, or c mirror planes 1696 ! onto the binary axis 1697 do isym=1,nsymconv 1698 if(t_axes(isym)==16)then 1699 if(abs(tnonsconv(direction,isym))>tol8)additional_info=2 1700 end if 1701 end do 1702 end if 1703 n_axest=(/0,0,0,0,1,0,0,1,1,0, 0,0,0,0,0,2,0,1,0,2, 0,0,0,0,0,0,0,0,0,0,0/) 1704 if(sum((n_axes-n_axest)**2)==0)then ! Spgroup 56 (Pccn) or 60 (Pbcn) 1705 write(std_out,*)' symspgr: 56 or 60' 1706 additional_info=1 1707 ! Select mirror plane n 1708 do isym=1,nsymconv 1709 if(t_axes(isym)==18)then 1710 ! Find direction of mirror plane 1711 if(symrelconv(1,1,isym)==-1)direction=1 1712 if(symrelconv(2,2,isym)==-1)direction=2 1713 if(symrelconv(3,3,isym)==-1)direction=3 1714 end if 1715 end do 1716 ! Examine the projection of the translation vector of the a, b, or c mirror planes 1717 ! onto the binary axis 1718 do isym=1,nsymconv 1719 if(t_axes(isym)==16)then 1720 if(abs(tnonsconv(direction,isym))<tol8)additional_info=2 1721 end if 1722 end do 1723 end if 1724 end if 1725 else if(brvltt==2)then 1726 ! In the few next lines, use additional_info as a flag 1727 additional_info=0 1728 ! If the bravais lattice is inner-centered 1729 if(nsymconv==8)then 1730 ! Test spgroup 23 (I222) or 24 (I2_{1}2_{1}2_{1}) 1731 n_axest=(/0,0,0,0,0,0,1,1,3,0, 0,0,0,0,0,0,0,0,0,3, 0,0,0,0,0,0,0,0,0,0,0/) 1732 if(sum((n_axes-n_axest)**2)==0) additional_info=1 1733 else if(nsymconv==24)then 1734 ! Test spgroup 197 (I23) or 199 (I2_{1}3) 1735 n_axest=(/0,0,0,0,0,0,1,1,3,16, 0,0,0,0,0,0,0,0,0,3, 0,0,0,0,0,0,0,0,0,0,0/) 1736 if(sum((n_axes-n_axest)**2)==0) additional_info=1 1737 end if 1738 if(additional_info==1)then 1739 write(std_out,*)' symspgr: (23 or 24) or (197 or 199)' 1740 ! Select the three binary axes (they might be 2 or 2_1 !) 1741 test_direction(:)=0 1742 do isym=1,nsymconv 1743 if(t_axes(isym)==20)then 1744 ! Find direction of axis 1745 do direction=1,3 1746 if(symrelconv(direction,direction,isym)==1)then 1747 test_direction(direction)=1 1748 if(abs(tnonsconv(direction,isym))<tol8)then 1749 vect(:,direction)=tnonsconv(:,isym) 1750 else 1751 vect(:,direction)=tnonsconv(:,isym)+half 1752 end if 1753 vect(:,direction)=vect(:,direction)-nint(vect(:,direction)-tol8) 1754 vect(direction,direction)=zero 1755 end if 1756 end do ! direction=1,3 1757 end if ! if binary axis 1758 end do ! isym 1759 if(test_direction(1)/=1 .or. test_direction(2)/=1 .and. test_direction(3)/=1)then 1760 write(msg, '(5a,3i4)' )& 1761 'For space groups 23, 24, 197 or 197, the three binary axes',ch10,& 1762 'are not equally partitioned along the x, y and z directions',ch10,& 1763 'test_direction(1:3)=',test_direction(:) 1764 ABI_BUG(msg) 1765 end if 1766 additional_info=1 1767 if(abs(vect(1,2)-vect(1,3))>tol8 .or. abs(vect(2,1)-vect(2,3))>tol8 .or. & 1768 abs(vect(3,1)-vect(3,2))>tol8) additional_info=2 1769 end if ! additional information are needed 1770 end if ! brvltt==1 1771 1772 if (brvltt==0 .or. brvltt==1) then ! Primitive 1773 call symlist_prim(additional_info,nsymconv,n_axes,spgroup) 1774 else if(brvltt==2)then 1775 call symlist_bcc(additional_info,nsymconv,n_axes,spgroup) 1776 else if(brvltt==3)then 1777 call symlist_fcc(nsymconv,n_axes,spgroup) 1778 else 1779 call symlist_others(brvltt,nsymconv,n_axes,spgroup) 1780 end if 1781 1782 if(spgroup==0) then 1783 write(msg, '(a,a,a,a,a)' )& 1784 'Could not find the space group.',ch10,& 1785 'This often happens when the user selects a restricted set of symmetries ',ch10,& 1786 'in the input file, instead of letting the code automatically find symmetries.' 1787 ABI_WARNING(msg) 1788 end if 1789 1790 spgorig=1 ; spgaxor=1 1791 call spgdata(brvsb,intsb,intsbl,ptintsb,ptschsb,schsb,spgaxor,spgroup,sporder,spgorig) 1792 1793 if(spgroup/=0)then 1794 write(msg, '(a,i4,2x,a,a,a,a,a)' ) ' symspgr: spgroup=',spgroup,trim(brvsb),trim(intsb),' (=',trim(schsb),')' 1795 call wrtout(std_out,msg) 1796 end if 1797 1798 if(bravais(1)==7)then 1799 write(msg, '(a)' ) ' symspgr: optical characteristics = isotropic ' 1800 call wrtout(std_out,msg) 1801 else if(bravais(1)==4 .or. bravais(1)==5 .or. bravais(1)==6)then 1802 write(msg, '(a)' ) ' symspgr: optical characteristics = uniaxial ' 1803 call wrtout(std_out,msg) 1804 ! Identify the first symmetry operation that is order 3, 4 or 6 1805 found=0 1806 do isym=1,nsym 1807 ! Proper rotations 1808 if( minval( abs( t_axes(isym)-(/10,12,14,22,23,24,25,26,27,28,29,30,31/) ))==0) then 1809 found=1 ; exit 1810 ! Improper symmetry operations 1811 else if( minval( abs( t_axes(isym)-(/1,2,3/) ))==0) then 1812 found=-1 ; exit 1813 end if 1814 end do 1815 if(found==-1 .or. found==1)then 1816 symrel_uni=symrel(:,:,isym) 1817 if(found==-1)symrel_uni=-symrel_uni 1818 ! Now, symrel_uni is a rotation of order 3, 4, 6, for which the axis must be identified 1819 ! It is actually the only eigenvector with eigenvalue 1. It can be found by cross products 1820 ! Subtract the unit matrix. 1821 do ii=1,3 1822 symrel_uni(ii,ii)=symrel_uni(ii,ii)-1 1823 end do 1824 found=0 1825 do ii=1,3 1826 jj=ii+1 ; if(jj==4)jj=1 1827 ! Cross product 1828 ivec1 = symrel_uni(ii,:); ivec2 = symrel_uni(jj,:) 1829 uniaxis = ivec1 .x. ivec2 1830 if(sum(uniaxis**2)/=0)then 1831 found=1 ; exit 1832 end if 1833 end do 1834 if(found==1)then 1835 ! Try to reduce the length, by an integer factor (try only primes 2, 3, 5, 7, 11) 1836 prime=(/2,3,5,7,11/) 1837 ii=1 1838 do while (ii<6) 1839 uniaxis_try=uniaxis/prime(ii) 1840 if(sum(abs(uniaxis_try*prime(ii)-uniaxis))==0)then 1841 uniaxis=uniaxis_try 1842 else 1843 ii=ii+1 1844 end if 1845 end do 1846 write(msg, '(a,3i4)' ) ' Optical axis (in reduced coordinates, real space ) :',uniaxis 1847 end if 1848 end if 1849 if(found==0)then 1850 write(msg, '(a)' ) ' However, the axis has not been found. Sorry for this.' 1851 end if 1852 call wrtout(std_out,msg) 1853 end if 1854 1855 ABI_FREE(determinant) 1856 ABI_FREE(shift) 1857 ABI_FREE(symrelconv) 1858 ABI_FREE(tnonsconv) 1859 ABI_FREE(t_axes) 1860 1861 DBG_EXIT("COLL") 1862 1863 end subroutine symspgr