TABLE OF CONTENTS
ABINIT/randomcellpos [ Functions ]
NAME
randomcellpos
FUNCTION
FIXME: This subroutine creates a unit cell with random atomic positions. It is
assumed that the cell parameters are given and fixed. Several methods are
used to generate the cell.
COPYRIGHT
Copyright (C) 2010-2017 ABINIT group (AHR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
natom=number of atoms npsp=number of pseudopotentials (needed for the dimension of znucl) ntypat=number of type of atoms random_atpos=input variable 0 no generation of random atomic potision 1 completely random atomic potisions 2 random atomic positions, avoiding too close atoms (prevent coming closer than a fraction of the sum of covalent radii) 3 same than 2 but also generates the rprim and acell randomly within some given ranges (angles between 50 and 130) ratsph(1:ntypat)=radius of the atomic sphere rprimd(3,3)=dimensional primitive translations in real space (bohr) typat(1:natom)= input variable giving the type of each atom znucl(1:npsp)=nuclear number of atom as specified in psp file
OUTPUT
xred(3,natom)=reduced dimensionless atomic coordinates
SIDE EFFECTS
NOTES
PARENTS
ingeo
CHILDREN
atomdata_from_znucl,mkrdim,xred2xcart
SOURCE
46 #if defined HAVE_CONFIG_H 47 #include "config.h" 48 #endif 49 50 #include "abi_common.h" 51 52 53 subroutine randomcellpos(natom,npsp,ntypat,random_atpos,ratsph,rprim,rprimd,typat,xred,znucl,acell) 54 55 use defs_basis 56 use m_profiling_abi 57 use m_errors 58 use m_atomdata 59 60 !This section has been created automatically by the script Abilint (TD). 61 !Do not modify the following lines by hand. 62 #undef ABI_FUNC 63 #define ABI_FUNC 'randomcellpos' 64 use interfaces_28_numeric_noabirule 65 use interfaces_41_geometry, except_this_one => randomcellpos 66 !End of the abilint section 67 68 implicit none 69 70 !Arguments ------------------------------------ 71 !scalars 72 integer,intent(in) :: natom,npsp,ntypat,random_atpos 73 !arrays 74 integer, intent(in) :: typat(natom) 75 real(dp),intent(in) :: ratsph(ntypat) 76 real(dp), intent(inout) :: rprim(3,3) 77 real(dp), intent(inout) :: rprimd(3,3) 78 real(dp), intent(inout) :: xred(3,natom) 79 real(dp), intent(in) :: znucl(npsp) 80 real(dp), intent(inout) :: acell(3) 81 82 !Local variables------------------------------- 83 integer :: iatom=0,ii,idum=-20 84 real(dp) :: rij(3), rijd(3), radiuscovi, radiuscovj, dist, rati, ratj, angdeg(3) 85 real(dp) :: cosang,aa,cc,a2 86 character(len=500) :: message 87 type(atomdata_t) :: atom 88 89 ! ************************************************************************* 90 91 !DEBUG 92 !For the time being, print rprimd to keep it as an argument, in spite of abirule checking. 93 !write (std_out,*) ' randomcellpos : enter' 94 !write(std_out,*)' rprimd=',rprimd 95 !write(std_out,*)' znucl=',znucl 96 !write(std_out,*)' typat=',typat 97 !write(std_out,*)' random_atpos=',random_atpos 98 !ENDDEBUG 99 100 if(random_atpos==2 .and. npsp/=ntypat)then 101 write(message, '(a,i5,2a,i5,a,i5,4a)' )& 102 & 'Input variable random_atpos= ',random_atpos,ch10,& 103 & 'However, the number of pseudopotentials ',npsp,', is not equal to the number of type of atoms ',ntypat,ch10,& 104 & 'The use of alchemical mixing cannot be combined with the constraint based on the mixing of covalent radii.',ch10,& 105 & 'Action : switch to another value of random_atpos.' 106 MSG_ERROR(message) 107 end if 108 109 !random_atpos = 0 Default value, no random initialisation 110 !random_atpos = 1 Fully random (Is it really useful ???) 111 !random_atpos = 2 Random, but the sum of the two covalent radii is 112 !less than the interatomic distance 113 !random_atpos = 3 Random, but the sum of the two (other type of) 114 !radii is less than the interatomic distance 115 !random_atpos = 4 Random, but the sum of the two pseudopotential 116 !radii is less than the interatomic distance 117 !random_atpos = 5 Random, but the interatomic distance must be bigger 118 !than the sum of 119 !some input variable (well, instead of defining a new variable, why 120 !not use ratsph ?) 121 !Right now we are not using a factor for the tested distance.. something to be done, after a new variable has been defined 122 123 if (random_atpos /= 0) then 124 select case (random_atpos) 125 case (1) 126 do ii=1,natom 127 xred(1,ii)=uniformrandom(idum) 128 xred(2,ii)=uniformrandom(idum) 129 xred(3,ii)=uniformrandom(idum) 130 end do 131 case (2) 132 iatom=0 133 do 134 iatom=iatom+1 135 xred(1,iatom)=uniformrandom(idum) 136 xred(2,iatom)=uniformrandom(idum) 137 xred(3,iatom)=uniformrandom(idum) 138 call atomdata_from_znucl(atom,znucl(typat(iatom))) 139 radiuscovi = atom%rcov 140 do ii=1,iatom-1 141 rij=xred(:,iatom)-xred(:,ii) 142 ! periodic boundary conditions 143 rij = rij - 0.5 144 rij = rij - anint (rij) 145 ! coming back to cube between (0,1) 146 rij = rij + 0.5 147 ! convert reduced coordinates to cartesian coordinates 148 call xred2xcart(1,rprimd,rijd,rij) 149 dist=dot_product(rijd,rijd) 150 call atomdata_from_znucl(atom,znucl(typat(ii))) 151 radiuscovj = atom%rcov 152 if (dist<(radiuscovj+radiuscovi)) then 153 iatom = iatom -1 154 EXIT 155 end if 156 end do 157 if (iatom>=natom) EXIT 158 end do 159 case(3) 160 iatom=0 161 do 162 iatom=iatom+1 163 xred(1,iatom)=uniformrandom(idum) 164 xred(2,iatom)=uniformrandom(idum) 165 xred(3,iatom)=uniformrandom(idum) 166 call atomdata_from_znucl(atom,znucl(typat(iatom))) 167 radiuscovi = atom%rcov 168 do ii=1,iatom-1 169 rij=xred(:,iatom)-xred(:,ii) 170 ! periodic boundary conditions 171 rij = rij - 0.5 172 rij = rij - anint (rij) 173 ! coming back to cube between (0,1) 174 rij = rij + 0.5 175 ! convert reduced coordinates to cartesian coordinates 176 call xred2xcart(1,rprimd,rijd,rij) 177 dist=dot_product(rijd,rijd) 178 call atomdata_from_znucl(atom,znucl(typat(ii))) 179 radiuscovj = atom%rcov 180 if (dist<(radiuscovj+radiuscovi)) then 181 iatom = iatom -1 182 EXIT 183 end if 184 end do 185 if (iatom>=natom) EXIT 186 end do 187 do ii=1,3 188 ! generates cells with angles between 60 and 120 degrees 189 angdeg(ii)=60_dp+uniformrandom(idum)*60.0_dp 190 end do 191 if (angdeg(1)+angdeg(2)+angdeg(3)>360._dp) then 192 angdeg(3)=360._dp-angdeg(1)-angdeg(2) 193 end if 194 ! check if angles are between the limits and create rprim 195 if( abs(angdeg(1)-angdeg(2))<tol12 .and. & 196 & abs(angdeg(2)-angdeg(3))<tol12 .and. & 197 & abs(angdeg(1)-90._dp)+abs(angdeg(2)-90._dp)+abs(angdeg(3)-90._dp)>tol12 )then 198 ! Treat the case of equal angles (except all right angles) : 199 ! generates trigonal symmetry wrt third axis 200 cosang=cos(pi*angdeg(1)/180.0_dp) 201 a2=2.0_dp/3.0_dp*(1.0_dp-cosang) 202 aa=sqrt(a2) 203 cc=sqrt(1.0_dp-a2) 204 rprim(1,1)=aa ; rprim(2,1)=0.0_dp ; rprim(3,1)=cc 205 rprim(1,2)=-0.5_dp*aa ; rprim(2,2)= sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,2)=cc 206 rprim(1,3)=-0.5_dp*aa ; rprim(2,3)=-sqrt(3.0_dp)*0.5_dp*aa ; rprim(3,3)=cc 207 ! DEBUG 208 ! write(std_out,*)' ingeo : angdeg=',angdeg(1:3) 209 ! write(std_out,*)' ingeo : aa,cc=',aa,cc 210 ! ENDDEBUG 211 else 212 ! Treat all the other cases 213 rprim(:,:)=0.0_dp 214 rprim(1,1)=1.0_dp 215 rprim(1,2)=cos(pi*angdeg(3)/180.0_dp) 216 rprim(2,2)=sin(pi*angdeg(3)/180.0_dp) 217 rprim(1,3)=cos(pi*angdeg(2)/180.0_dp) 218 rprim(2,3)=(cos(pi*angdeg(1)/180.0_dp)-rprim(1,2)*rprim(1,3))/rprim(2,2) 219 rprim(3,3)=sqrt(1.0_dp-rprim(1,3)**2-rprim(2,3)**2) 220 end if 221 ! generate acell 222 aa=zero 223 do ii=1,npsp 224 aa=znucl(ii) 225 end do 226 do ii=1,3 227 acell(ii)=aa+uniformrandom(idum)*4.0 228 end do 229 call mkrdim(acell,rprim,rprimd) 230 case(4) 231 write(std_out,*) 'Not implemented yet' 232 case(5) 233 iatom=0 234 do 235 iatom=iatom+1 236 xred(1,iatom)=uniformrandom(idum) 237 xred(2,iatom)=uniformrandom(idum) 238 xred(3,iatom)=uniformrandom(idum) 239 rati=ratsph(typat(iatom)) 240 do ii=1,iatom-1 241 ratj=ratsph(typat(ii)) 242 ! apply periodic boundary conditions 243 rij=(xred(:,iatom)-xred(:,ii))-0.5 244 rij = rij - ANINT ( rij ) 245 rij = rij + 0.5 246 call xred2xcart(natom,rprimd,rijd,rij) 247 dist=dot_product(rijd,rijd) 248 if (dist<(rati+ratj)) EXIT 249 end do 250 if (iatom==natom) EXIT 251 if (ii<(iatom-1)) iatom=iatom-1 252 end do 253 end select 254 end if 255 256 end subroutine randomcellpos