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