TABLE OF CONTENTS
ABINIT/bonds_lgth_angles [ Functions ]
NAME
bonds_lgth_angles
FUNCTION
From list of coordinates and primitive translations, output a list of bonds lengths and bond angles.
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (DCA, 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
coordn = maximum coordination number to be taken into account fnameabo_app_geo=name of file for _GEO data natom = number of atoms in unit cell ntypat = number of types of atoms in unit cell. rprimd(3,3) = real space dimensional primitive translations (bohr) typat(natom) = type integer for each atom in cell znucl(ntypat)= real(dp), atomic number of atom type xred(3,natom)= reduced coordinates of atoms
OUTPUT
data written in file fnameabo_app_geo
NOTES
The tolerance tol8 aims at giving a machine-independent ordering. (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)
PARENTS
outscfcv
CHILDREN
atomdata_from_znucl,wrtout,xred2xcart
SOURCE
42 #if defined HAVE_CONFIG_H 43 #include "config.h" 44 #endif 45 46 #include "abi_common.h" 47 48 subroutine bonds_lgth_angles(coordn,fnameabo_app_geo,natom,ntypat,rprimd,typat,xred,znucl) 49 50 use defs_basis 51 use m_errors 52 use m_profiling_abi 53 use m_atomdata 54 use m_errors 55 56 use m_io_tools, only : open_file 57 58 !This section has been created automatically by the script Abilint (TD). 59 !Do not modify the following lines by hand. 60 #undef ABI_FUNC 61 #define ABI_FUNC 'bonds_lgth_angles' 62 use interfaces_14_hidewrite 63 use interfaces_41_geometry, except_this_one => bonds_lgth_angles 64 !End of the abilint section 65 66 implicit none 67 68 !Arguments ------------------------------------ 69 !scalars 70 integer,intent(in) :: coordn,natom,ntypat 71 character(len=*),intent(in) :: fnameabo_app_geo 72 !arrays 73 integer,intent(in) :: typat(natom) 74 real(dp),intent(in) :: rprimd(3,3),znucl(ntypat) 75 real(dp),intent(inout) :: xred(3,natom) 76 77 !Local variables------------------------------- 78 !scalars 79 integer :: done,ia,ib,ic,ii,ineighb,jneighb,mneighb,mu,ndig,nu,t1,t2,t3,tmax,temp_unit 80 real(dp) :: adotb,asq,bsq,co,length,sq,thdeg 81 !real(dp)u1,u2,u3,v1,v2,v3 82 character(len=500) :: message 83 type(atomdata_t) :: atom 84 !arrays 85 integer,allocatable :: list_neighb(:,:,:) 86 real(dp) :: bab(3),bac(3),dif(3),rmet(3,3) 87 real(dp),allocatable :: sqrlength(:),xangst(:,:),xcart(:,:) 88 character(len=8),allocatable :: iden(:) 89 90 ! ************************************************************************* 91 92 !Initialize the file 93 write(message, '(a,a,a)' )' bonds_lgth_angles : about to open file ',trim(fnameabo_app_geo),ch10 94 call wrtout(std_out,message,'COLL'); call wrtout(ab_out,message,'COLL') 95 96 if (open_file(fnameabo_app_geo,message,newunit=temp_unit,status='unknown',form='formatted') /= 0) then 97 MSG_ERROR(message) 98 end if 99 rewind(temp_unit) 100 101 write(message, '(a,a)' ) ch10,' ABINIT package : GEO file ' 102 call wrtout(temp_unit,message,'COLL') 103 104 !Compute maximum number of neighbors is the neighbor list, 105 !from the indicative coordination number 106 !Note : the following formula includes next nearest neighbors, but not others 107 mneighb=1+coordn+coordn*(coordn-1) 108 109 write(message, '(a,a,i2,a,a,i4,a,a,a,i4,a)' ) ch10,& 110 & ' Maximal coordination number, as estimated by the user : ',coordn,ch10,& 111 & ' giving a maximum of ',coordn*coordn,& 112 & ' nearest neighbors and next nearest neighbors, ',ch10,& 113 & ' and ',(coordn*(coordn-1))/2,& 114 & ' distinct angles between nearest neighbors' 115 call wrtout(temp_unit,message,'COLL') 116 117 !Compute metric tensor in real space rmet 118 do nu=1,3 119 do mu=1,3 120 rmet(mu,nu)=rprimd(1,mu)*rprimd(1,nu)+& 121 & rprimd(2,mu)*rprimd(2,nu)+& 122 & rprimd(3,mu)*rprimd(3,nu) 123 end do 124 end do 125 126 write(message, '(a,a)' )ch10,' Primitive vectors of the periodic cell (bohr)' 127 call wrtout(temp_unit,message,'COLL') 128 do nu=1,3 129 write(message, '(1x,a,i1,a,3f10.5)' ) ' R(',nu,')=',rprimd(:,nu) 130 call wrtout(temp_unit,message,'COLL') 131 end do 132 133 write(message, '(a,a)' ) ch10,& 134 & ' Atom list Reduced coordinates Cartesian coordinates (bohr)' 135 call wrtout(temp_unit,message,'COLL') 136 137 !Set up a list of character identifiers for all atoms : iden(ia) 138 ABI_ALLOCATE(iden,(natom)) 139 iden(:)=' ' 140 do ia=1,natom 141 ndig=int(log10(dble(ia)+0.5d0))+1 142 call atomdata_from_znucl(atom,znucl(typat(ia))) 143 if(ndig==1) write(iden(ia), '(a,a,i1,a)' ) atom%symbol,'(',ia,') ' 144 if(ndig==2) write(iden(ia), '(a,a,i2,a)' ) atom%symbol,'(',ia,') ' 145 if(ndig==3) write(iden(ia), '(a,a,i3,a)' ) atom%symbol,'(',ia,') ' 146 if(ndig==4) write(iden(ia), '(a,a,i4,a)' ) atom%symbol,'(',ia,')' 147 if(ndig>4)then 148 close(temp_unit) 149 write(message, '(a,i8,a,a)' )& 150 & 'bonds_lgth_angles cannot handle more than 9999 atoms, while natom=',natom,ch10,& 151 & 'Action : decrease natom, or contact ABINIT group.' 152 MSG_BUG(message) 153 end if 154 end do 155 156 !Compute cartesian coordinates, and print reduced and cartesian coordinates 157 !then print coordinates in angstrom, with the format neede for xmol 158 ABI_ALLOCATE(xangst,(3,natom)) 159 ABI_ALLOCATE(xcart,(3,natom)) 160 call xred2xcart(natom,rprimd,xcart,xred) 161 xangst(:,:)=xcart(:,:)*Bohr_Ang 162 163 do ia=1,natom 164 write(message, '(a,a,3f10.5,a,3f10.5)' ) & 165 & ' ',iden(ia),(xred(ii,ia)+tol10,ii=1,3),& 166 & ' ',(xcart(ii,ia)+tol10,ii=1,3) 167 call wrtout(temp_unit,message,'COLL') 168 end do 169 170 write(message, '(a,a,a,a,i4,a)' )ch10,& 171 & ' XMOL data : natom, followed by cartesian coordinates in Angstrom',& 172 & ch10,ch10,natom,ch10 173 call wrtout(temp_unit,message,'COLL') 174 175 do ia=1,natom 176 call atomdata_from_znucl(atom,znucl(typat(ia))) 177 write(message, '(a,a,3f10.5)' )' ',atom%symbol,xangst(1:3,ia) 178 call wrtout(temp_unit,message,'COLL') 179 end do 180 181 ABI_DEALLOCATE(xangst) 182 ABI_DEALLOCATE(xcart) 183 184 ABI_ALLOCATE(list_neighb,(0:mneighb+1,4,2)) 185 ABI_ALLOCATE(sqrlength,(0:mneighb+1)) 186 187 !Compute list of neighbors 188 do ia=1,natom 189 190 write(message, '(a,a,a,a,a,a,a,a,a)' ) ch10,'===========',& 191 & '=====================================================================',& 192 & ch10,' ',iden(ia),ch10,ch10,' Bond lengths ' 193 call wrtout(temp_unit,message,'COLL') 194 195 ! Search other atoms for bonds, but must proceed 196 ! in such a way to consider a search box sufficiently large, 197 ! so increase the size of the search box until the 198 ! final bond length list do not change 199 do tmax=0,5 200 201 ! Set initial list of neighbors to zero, 202 ! and initial square of bond lengths to a very large number. 203 ! Note that the dimension is larger than neighb to ease 204 ! the later sorting : neighbors 0 and neighb+1 are non-existent, while 205 ! neighbor 1 will be the atom itself ... 206 list_neighb(0:mneighb+1,1:4,1)=0 207 sqrlength(1:mneighb+1)=huge(0.0d0) 208 sqrlength(0)=-1.0d0 209 210 ! Here search on all atoms inside the box defined by tmax 211 do ib=1,natom 212 do t3=-tmax,tmax 213 do t2=-tmax,tmax 214 do t1=-tmax,tmax 215 dif(1)=xred(1,ia)-(xred(1,ib)+dble(t1)) 216 dif(2)=xred(2,ia)-(xred(2,ib)+dble(t2)) 217 dif(3)=xred(3,ia)-(xred(3,ib)+dble(t3)) 218 sq=rsdot(dif(1),dif(2),dif(3),dif(1),dif(2),dif(3),rmet) 219 220 ! Insert the atom at the proper place in the neighbor list. 221 do ineighb=mneighb,0,-1 222 ! Note the tolerance 223 if(sq+tol8>sqrlength(ineighb))then 224 sqrlength(ineighb+1)=sq 225 list_neighb(ineighb+1,1,1)=ib 226 list_neighb(ineighb+1,2,1)=t1 227 list_neighb(ineighb+1,3,1)=t2 228 list_neighb(ineighb+1,4,1)=t3 229 ! DEBUG 230 ! if(ineighb/=mneighb)then 231 ! write(std_out,*)' ' 232 ! do ii=1,mneighb 233 ! write(std_out,*)ii,sqrlength(ii) 234 ! end do 235 ! end if 236 ! ENDDEBUG 237 exit 238 else 239 sqrlength(ineighb+1)=sqrlength(ineighb) 240 list_neighb(ineighb+1,1:4,1)=list_neighb(ineighb,1:4,1) 241 end if 242 end do 243 244 end do 245 end do 246 end do 247 ! end ib loop: 248 end do 249 250 ! Now, check that the box defined by tmax was large enough : 251 ! require the present and old lists to be the same 252 done=0 253 254 if(tmax>0)then 255 done=1 256 do ineighb=1,mneighb 257 ! DEBUG 258 ! write(std_out,'(5i5,f12.5)' )ineighb,list_neighb(ineighb,1:4,1),& 259 ! & sqrlength(ineighb) 260 ! write(std_out,'(5i5)' )ineighb,list_neighb(ineighb,1:4,2) 261 ! ENDDEBUG 262 if( list_neighb(ineighb,1,1)/=list_neighb(ineighb,1,2) .or. & 263 & list_neighb(ineighb,2,1)/=list_neighb(ineighb,2,2) .or. & 264 & list_neighb(ineighb,3,1)/=list_neighb(ineighb,3,2) .or. & 265 & list_neighb(ineighb,4,1)/=list_neighb(ineighb,4,2) )then 266 done=0 267 end if 268 end do 269 end if 270 271 ! If done==1, then one can exit the loop : the correct list of 272 ! neighbors is contained in list_neighb(1:neighb,1:4,1), 273 ! with the first neighbor being the atom itself 274 if(done==1)exit 275 276 ! If the work is not done, while tmax==5, then there is a problem . 277 if(tmax==5)then 278 close(temp_unit) 279 write(message, '(2a)' )& 280 & 'Did not succeed to generate a reliable list of bonds ',& 281 & 'since tmax is exceeded.' 282 MSG_BUG(message) 283 end if 284 285 ! Copy the new list into the old list. 286 list_neighb(1:mneighb,1:4,2)=list_neighb(1:mneighb,1:4,1) 287 288 ! Loop on tmax (note that there are exit instruction inside the loop) 289 end do 290 291 292 293 ! Output the bond list 294 do ineighb=2,mneighb 295 ib=list_neighb(ineighb,1,1) 296 length=sqrt(sqrlength(ineighb)) 297 write(message, '(a,a,a,a,3i2,t27,a,f10.5,a,f9.5,a)' )& 298 & ' ',trim(iden(ia)),' - ',trim(iden(ib)),& 299 & list_neighb(ineighb,2:4,1),'bond length is ',& 300 & length,' bohr ( or ',Bohr_Ang*length,' Angst.)' 301 call wrtout(temp_unit,message,'COLL') 302 end do 303 304 ! Output the angle list 305 if(coordn>1)then 306 307 write(message, '(a,a)' ) ch10,' Bond angles ' 308 call wrtout(temp_unit,message,'COLL') 309 310 do ineighb=2,coordn 311 do jneighb=ineighb+1,coordn+1 312 313 ib=list_neighb(ineighb,1,1) 314 ic=list_neighb(jneighb,1,1) 315 do mu=1,3 316 bab(mu)=xred(mu,ib)+dble(list_neighb(ineighb,1+mu,1))-xred(mu,ia) 317 bac(mu)=xred(mu,ic)+dble(list_neighb(jneighb,1+mu,1))-xred(mu,ia) 318 end do 319 asq=rsdot(bab(1),bab(2),bab(3),bab(1),bab(2),bab(3),rmet) 320 bsq=rsdot(bac(1),bac(2),bac(3),bac(1),bac(2),bac(3),rmet) 321 adotb=rsdot(bab(1),bab(2),bab(3),bac(1),bac(2),bac(3),rmet) 322 co=adotb/sqrt(asq*bsq) 323 if( abs(co)-1.0d0 >= 0.0d0 )then 324 if( abs(co)-1.0d0 <= 1.0d-12 )then 325 ! Allows for a small numerical inaccuracy 326 thdeg=0.0d0 327 if(co < 0.0d0) thdeg=180.0d0 328 else 329 MSG_BUG('the evaluation of the angle is wrong.') 330 end if 331 else 332 thdeg=acos(co)*180.d0*piinv 333 end if 334 335 write(message, '(a,a,3i2,a,a,a,a,3i2,t44,a,f13.5,a)' )& 336 & ' ',trim(iden(ib)),list_neighb(ineighb,2:4,1),' - ',& 337 & trim(iden(ia)),' - ',trim(iden(ic)),& 338 & list_neighb(jneighb,2:4,1),'bond angle is ',thdeg,' degrees ' 339 call wrtout(temp_unit,message,'COLL') 340 end do 341 end do 342 343 end if 344 end do ! End big ia loop: 345 346 ABI_DEALLOCATE(iden) 347 ABI_DEALLOCATE(list_neighb) 348 ABI_DEALLOCATE(sqrlength) 349 350 close(temp_unit) 351 352 contains 353 354 function rsdot(u1,u2,u3,v1,v2,v3,rmet) 355 356 357 !This section has been created automatically by the script Abilint (TD). 358 !Do not modify the following lines by hand. 359 #undef ABI_FUNC 360 #define ABI_FUNC 'rsdot' 361 !End of the abilint section 362 363 real(dp) :: rsdot 364 real(dp),intent(in) :: u1,u2,u3,v1,v2,v3 365 real(dp),intent(in) :: rmet(3,3) 366 rsdot=rmet(1,1)*u1*v1+rmet(2,1)*u2*v1+& 367 & rmet(3,1)*u3*v1+rmet(1,2)*u1*v2+rmet(2,2)*u2*v2+& 368 & rmet(3,2)*u3*v2+rmet(1,3)*u1*v3+rmet(2,3)*u2*v3+rmet(3,3)*u3*v3 369 end function rsdot 370 371 end subroutine bonds_lgth_angles