TABLE OF CONTENTS
ABINIT/smallprim [ Functions ]
NAME
smallprim
FUNCTION
Find the smallest possible primitive vectors for an input lattice This algorithm is not as restrictive as the conditions mentioned at p.740 of the international tables for crystallography (1983). The final vectors form a right-handed basis, while their sign and ordering is chosen such as to maximize the overlap with the original vectors in order.
COPYRIGHT
Copyright (C) 2000-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
rprimd(3,3)=primitive vectors
OUTPUT
metmin(3,3)=metric for the new (minimal) primitive vectors minim(3,3)=minimal primitive translations
NOTES
The routine might as well be defined without metmin as argument, but it is more convenient to have it
PARENTS
getkgrid,symlatt,testkgrid
CHILDREN
metric
SOURCE
41 #if defined HAVE_CONFIG_H 42 #include "config.h" 43 #endif 44 45 #include "abi_common.h" 46 47 48 subroutine smallprim(metmin,minim,rprimd) 49 50 use defs_basis 51 use m_errors 52 use m_profiling_abi 53 54 use m_geometry, only : metric 55 56 !This section has been created automatically by the script Abilint (TD). 57 !Do not modify the following lines by hand. 58 #undef ABI_FUNC 59 #define ABI_FUNC 'smallprim' 60 !End of the abilint section 61 62 implicit none 63 64 !Arguments ------------------------------------ 65 !arrays 66 real(dp),intent(in) :: rprimd(3,3) 67 real(dp),intent(out) :: metmin(3,3),minim(3,3) 68 69 !Local variables------------------------------- 70 !scalars 71 integer :: ia,ib,ii,itrial,minimal 72 integer :: iiter, maxiter = 100000 73 real(dp) :: determinant,length2,metsum,ucvol 74 character(len=500) :: message 75 !arrays 76 integer :: nvecta(3),nvectb(3) 77 real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3),scprod(3),tmpvect(3) 78 79 !************************************************************************** 80 81 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 82 83 nvecta(1)=2 ; nvectb(1)=3 84 nvecta(2)=1 ; nvectb(2)=3 85 nvecta(3)=1 ; nvectb(3)=2 86 87 minim(:,:)=rprimd(:,:) 88 metmin(:,:)=rmet(:,:) 89 90 !DEBUG 91 !write(std_out,*)' smallprim : starting values, rprim ' 92 !write(std_out,'(3f16.8)' )rprimd(:,1) 93 !write(std_out,'(3f16.8)' )rprimd(:,2) 94 !write(std_out,'(3f16.8)' )rprimd(:,3) 95 !write(std_out,*)' smallprim : starting values, rmet ' 96 !write(std_out,'(3f16.8)' )rmet(:,1) 97 !write(std_out,'(3f16.8)' )rmet(:,2) 98 !write(std_out,'(3f16.8)' )rmet(:,3) 99 !ENDDEBUG 100 101 !Note this loop without index 102 do iiter = 1, maxiter 103 104 ! Will exit if minimal=1 is still valid after a trial 105 ! to reduce the vectors of each of the three pairs 106 minimal=1 107 108 do itrial=1,3 109 110 ia=nvecta(itrial) ; ib=nvectb(itrial) 111 ! Make sure the scalar product is negative 112 if(metmin(ia,ib)>tol8)then 113 minim(:,ia)=-minim(:,ia) 114 metmin(ia,ib)=-metmin(ia,ib) ; metmin(ib,ia)=-metmin(ib,ia) 115 metmin(ia,itrial)=-metmin(ia,itrial) 116 metmin(itrial,ia)=-metmin(itrial,ia) 117 end if 118 ! Compute the length of the sum vector 119 length2=metmin(ia,ia)+2*metmin(ia,ib)+metmin(ib,ib) 120 ! Replace the first vector by the sum vector if the latter is smaller 121 if(length2/metmin(ia,ia) < one-tol8)then 122 minim(:,ia)=minim(:,ia)+minim(:,ib) 123 metmin(ia,ia)=length2 124 metmin(ia,ib)=metmin(ia,ib)+metmin(ib,ib) 125 metmin(ia,itrial)=metmin(ia,itrial)+metmin(ib,itrial) 126 metmin(ib,ia)=metmin(ia,ib) 127 metmin(itrial,ia)=metmin(ia,itrial) 128 minimal=0 129 ! Replace the second vector by the sum vector if the latter is smaller 130 else if(length2/metmin(ib,ib) < one-tol8)then 131 minim(:,ib)=minim(:,ia)+minim(:,ib) 132 metmin(ib,ib)=length2 133 metmin(ia,ib)=metmin(ia,ib)+metmin(ia,ia) 134 metmin(itrial,ib)=metmin(itrial,ib)+metmin(itrial,ia) 135 metmin(ib,ia)=metmin(ia,ib) 136 metmin(ib,itrial)=metmin(itrial,ib) 137 minimal=0 138 end if 139 140 end do 141 142 if(minimal==1)exit 143 144 end do 145 146 if (iiter >= maxiter) then 147 write(message,'(a,i0,a)') & 148 & 'the loop has failed to find a set of minimal vectors in ',maxiter,' iterations.' 149 MSG_BUG(message) 150 end if 151 152 !At this stage, the three vectors have angles between each other that are 153 !comprised between 90 and 120 degrees. It might still be that minus the vector 154 !that is the sum of the three vectors is smaller than the longest of these vectors 155 do iiter = 1, maxiter 156 157 ! Will exit if minimal=1 is still valid after a trial 158 ! to replace each of the three vectors by minus the summ of the three vectors 159 minimal=1 160 metsum=sum(metmin(:,:)) 161 do itrial=1,3 162 ia=nvecta(itrial) ; ib=nvectb(itrial) 163 if(metmin(ia,ia)/metsum > one + tol8)then 164 minim(:,ia)=-minim(:,1)-minim(:,2)-minim(:,3) 165 metmin(ia,ib)=-sum(metmin(:,ib)) 166 metmin(ia,itrial)=-sum(metmin(:,itrial)) 167 metmin(ia,ia)=metsum 168 metmin(ib,ia)=metmin(ia,ib) 169 metmin(itrial,ia)=metmin(ia,itrial) 170 minimal=0 171 end if 172 end do 173 174 if(minimal==1)exit 175 176 end do 177 178 if (iiter >= maxiter) then 179 write(message, '(a,i0,a)') & 180 & 'the second loop has failed to find a set of minimal vectors in ',maxiter, 'iterations.' 181 MSG_BUG(message) 182 end if 183 184 !DEBUG 185 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',& 186 !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3) 187 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' minim =',& 188 !& minim(:,1),ch10,minim(:,2),ch10,minim(:,3) 189 !ENDDEBUG 190 191 !DEBUG 192 !Change sign of the third vector if not right-handed basis 193 !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+& 194 !& minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+& 195 !& minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3)) 196 !write(std_out,*)' smallprim: determinant=',determinant 197 !ENDDEBUG 198 199 !Choose the first vector 200 !Compute the scalar product of the three minimal vectors 201 !with the first original vector 202 scprod(:)=zero 203 do ii=1,3 204 scprod(:)=scprod(:)+minim(ii,:)*rprimd(ii,1) 205 end do 206 !Determine the vector with the maximal absolute overlap 207 itrial=1 208 if(abs(scprod(2))>abs(scprod(1))+tol8)itrial=2 209 if(abs(scprod(3))>abs(scprod(itrial))+tol8)itrial=3 210 !Switch the vectors if needed 211 if(itrial/=1)then 212 tmpvect(:)=minim(:,1) 213 minim(:,1)=minim(:,itrial) 214 minim(:,itrial)=tmpvect(:) 215 end if 216 !Choose the sign 217 if(scprod(itrial)<tol8)minim(:,1)=-minim(:,1) 218 219 !DEBUG 220 !Change sign of the third vector if not right-handed basis 221 !determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+& 222 !& minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+& 223 !& minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3)) 224 !write(std_out,*)' smallprim: determinant=',determinant 225 !ENDDEBUG 226 227 !Choose the second vector 228 !Compute the scalar product of the second and third minimal vectors 229 !with the second original vector 230 scprod(2:3)=zero 231 do ii=1,3 232 scprod(2:3)=scprod(2:3)+minim(ii,2:3)*rprimd(ii,2) 233 end do 234 !Determine the vector with the maximal absolute overlap 235 itrial=2 236 if(abs(scprod(3))>abs(scprod(2))+tol8)itrial=3 237 !Switch the vectors if needed 238 if(itrial/=2)then 239 tmpvect(:)=minim(:,2) 240 minim(:,2)=minim(:,itrial) 241 minim(:,itrial)=tmpvect(:) 242 end if 243 !Choose the sign 244 if(scprod(itrial)<tol8)minim(:,2)=-minim(:,2) 245 246 !Change sign of the third vector if not right-handed basis 247 determinant=minim(1,1)*(minim(2,2)*minim(3,3)-minim(3,2)*minim(2,3))+& 248 & minim(2,1)*(minim(3,2)*minim(1,3)-minim(1,2)*minim(3,3))+& 249 & minim(3,1)*(minim(1,2)*minim(2,3)-minim(2,2)*minim(1,3)) 250 if(determinant<-tol8)minim(:,3)=-minim(:,3) 251 if(abs(determinant)<tol8)then 252 MSG_BUG('minim gives vanishing unit cell volume.') 253 end if 254 255 !Final computation of metmin 256 do ii=1,3 257 metmin(ii,:)=minim(1,ii)*minim(1,:)+& 258 & minim(2,ii)*minim(2,:)+& 259 & minim(3,ii)*minim(3,:) 260 end do 261 262 !DEBUG 263 !write(std_out,'(a,3es14.6,a,3es14.6,a,3es14.6)')' rprimd=',& 264 !& rprimd(:,1),ch10,rprimd(:,2),ch10,rprimd(:,3) 265 !write(std_out,'(a,3es16.8,a,3es16.8,a,3es16.8)')' minim =',& 266 !& minim(:,1),ch10,minim(:,2),ch10,minim(:,3) 267 !write(std_out,'(a,3es16.8,a,3es16.8,a,3es16.8)')' metmin =',& 268 !& metmin(:,1),ch10,metmin(:,2),ch10,metmin(:,3) 269 !ENDDEBUG 270 271 end subroutine smallprim