TABLE OF CONTENTS
- ABINIT/m_polynomial_term
- m_polynomial_term/polynomial_term_free
- m_polynomial_term/polynomial_term_init
- m_polynomial_term/polynomial_term_type
- m_polynomial_term/terms_compare
ABINIT/m_polynomial_term [ Functions ]
NAME
m_polynomial_term
FUNCTION
Module with the datatype polynomial terms
COPYRIGHT
Copyright (C) 2010-2024 ABINIT group (AM) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
SOURCE
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_polynomial_term 26 27 use defs_basis 28 use m_errors 29 use m_abicore 30 31 implicit none 32 33 public :: polynomial_term_init 34 public :: polynomial_term_free
m_polynomial_term/polynomial_term_free [ Functions ]
[ Top ] [ m_polynomial_term ] [ Functions ]
NAME
polynomial_term_free
FUNCTION
Free polynomial_term
INPUTS
polynomial_term<type(polynomial_term)> = datatype to free
OUTPUT
polynomial_term<type(polynomial_term)> = datatype to free
SOURCE
308 subroutine polynomial_term_free(polynomial_term) 309 310 implicit none 311 312 !Arguments ------------------------------------ 313 !scalars 314 !arrays 315 type(polynomial_term_type), intent(inout) :: polynomial_term 316 !Local variables------------------------------- 317 !scalar 318 !arrays 319 320 ! ************************************************************************* 321 322 polynomial_term%ndisp = 0 323 polynomial_term%nstrain = 0 324 polynomial_term%weight = zero 325 326 if(allocated(polynomial_term%atindx))then 327 polynomial_term%atindx(:,:) = 0 328 ABI_FREE(polynomial_term%atindx) 329 end if 330 331 if(allocated(polynomial_term%cell))then 332 polynomial_term%cell(:,:,:) = 0 333 ABI_FREE(polynomial_term%cell) 334 end if 335 336 if(allocated(polynomial_term%direction))then 337 polynomial_term%direction(:) = 0 338 ABI_FREE(polynomial_term%direction) 339 end if 340 341 if(allocated(polynomial_term%power_disp))then 342 polynomial_term%power_disp(:) = 0 343 ABI_FREE(polynomial_term%power_disp) 344 end if 345 346 if(allocated(polynomial_term%power_strain))then 347 polynomial_term%power_strain(:) = 0 348 ABI_FREE(polynomial_term%power_strain) 349 end if 350 351 if(allocated(polynomial_term%strain))then 352 polynomial_term%strain(:) = 0 353 ABI_FREE(polynomial_term%strain) 354 end if 355 356 end subroutine polynomial_term_free
m_polynomial_term/polynomial_term_init [ Functions ]
[ Top ] [ m_polynomial_term ] [ Functions ]
NAME
polynomial_term_init
FUNCTION
Initialize a polynomial_term_init for given set of displacements/strain
INPUTS
atindx(2) = Indexes of the atoms a and b in the unit cell
cell(3,2) = Indexes of the cell of the atom a and b
direction = direction of the perturbation => 1,2,3 for atomic displacement
ndisp = Number of displacement for this terms
nstrain = Number of strain for this terms
power_disp = power_disp of the displacement 2 (X_z-O_z)^2 or 1 for (X_y-O_y)^1
power_strain = power_strain of the strain 2 (\eta)^2 or 1 for (\eta)^1
strain(nstrain) = index of strain, 1 2 3 4 5 or 6
weight = Weight of the term
check = optional,logical => if TRUE, the term will be check, same displacement/strain
are gathered in the same displacement but with an higher
power_disp. For example:
((Sr_y-O1_y)^1(Sr_y-O1_y)^1 => (Sr_y-O1_y)^2)
if FALSE, default, do nothing
OUTPUT
polynomial_term<type(polynomial_term)> = polynomial_term datatype is now initialized
SOURCE
122 subroutine polynomial_term_init(atindx,cell,direction,ndisp,nstrain,polynomial_term,power_disp,& 123 & power_strain,strain,weight,check) 124 125 implicit none 126 127 !Arguments ------------------------------------ 128 !scalars 129 integer, intent(in) :: ndisp,nstrain 130 real(dp),intent(in) :: weight 131 logical,optional,intent(in) :: check 132 !arrays 133 integer, intent(in) :: atindx(2,ndisp) 134 integer, intent(in) :: cell(3,2,ndisp) 135 integer, intent(in) :: direction(ndisp),power_disp(ndisp) 136 integer, intent(in) :: strain(nstrain),power_strain(nstrain) 137 type(polynomial_term_type), intent(out) :: polynomial_term 138 !Local variables------------------------------- 139 !scalar 140 integer :: idisp1,idisp2,ndisp_tmp,nstrain_tmp 141 logical :: check_in 142 !arrays 143 integer :: power_disp_tmp(ndisp),power_strain_tmp(nstrain) 144 character(500) :: msg 145 146 ! ************************************************************************* 147 148 !Do some checks 149 if (size(atindx,2) /= ndisp) then 150 write(msg,'(a)')' atindx and ndisp have not the same size' 151 ABI_ERROR(msg) 152 end if 153 154 if (size(cell,3) /= ndisp) then 155 write(msg,'(a)')' cell and ndisp have not the same size' 156 ABI_ERROR(msg) 157 end if 158 159 if (size(direction) /= ndisp) then 160 write(msg,'(a)')' direction and ndisp have not the same size' 161 ABI_ERROR(msg) 162 end if 163 164 if (size(power_disp) /= ndisp) then 165 write(msg,'(a)')' power_disp and ndisp have not the same size' 166 ABI_ERROR(msg) 167 end if 168 169 if (size(power_strain) /= nstrain) then 170 write(msg,'(a)')' power_strain and nstrain have not the same size' 171 ABI_ERROR(msg) 172 end if 173 174 if (size(strain) /= nstrain) then 175 write(msg,'(a)')' strain and nstrain have not the same size' 176 ABI_ERROR(msg) 177 end if 178 179 ! FIXME: hexu: check why this does not work? 180 ! if (ndisp>1) then 181 ! print *, "atinx(1, 1:ndisp)", atindx(1, 1:ndisp) 182 ! if (.not. all(atindx(1, 1:ndisp)-atindx(1,1)==0)) then 183 ! write(msg,'(a)')' Not all displacement pairs start with the same atom.' 184 ! print *, "atinx(1, 1:ndisp)", atindx(1, 1:ndisp) 185 ! ABI_ERROR(msg) 186 ! end if 187 ! end if 188 189 !First free datatype before init 190 call polynomial_term_free(polynomial_term) 191 check_in = .false. 192 193 !Copy the power array before check 194 power_disp_tmp(:) = power_disp(:) 195 power_strain_tmp(:) = power_strain(:) 196 197 if(present(check)) check_in = check 198 199 if(check_in)then 200 !Check if displacement are identical, in this case 201 !increase the power_disp 202 do idisp1=1,ndisp 203 do idisp2=idisp1,ndisp 204 if (idisp1/=idisp2.and.& 205 & atindx(1,idisp1) == atindx(1,idisp2).and.& 206 & atindx(2,idisp1) == atindx(2,idisp2).and.& 207 & direction(idisp1) == direction(idisp2).and.& 208 & all(cell(:,1,idisp1)==cell(:,1,idisp2)).and.& 209 & all(cell(:,2,idisp1)==cell(:,2,idisp2)).and.& 210 & power_disp_tmp(idisp2) > 0 )then 211 power_disp_tmp(idisp1) = power_disp_tmp(idisp1) + 1 212 power_disp_tmp(idisp2) = 0 213 end if 214 end do 215 end do 216 217 ! Count the number of power_disp avec the previous check 218 ! or just remove the power_disp equal to zero 219 ndisp_tmp = 0 220 do idisp1=1,ndisp 221 if(power_disp_tmp(idisp1) > zero) then 222 ndisp_tmp = ndisp_tmp + 1 223 end if 224 end do 225 226 !Check if strain are identical, in this case 227 !increase the power_strain 228 do idisp1=1,nstrain 229 do idisp2=idisp1,nstrain 230 if (idisp1/=idisp2.and.& 231 & strain(idisp1) == strain(idisp2).and.& 232 & power_strain_tmp(idisp2) > 0 )then 233 power_strain_tmp(idisp1) = power_strain_tmp(idisp1) + 1 234 power_strain_tmp(idisp2) = 0 235 end if 236 end do 237 end do 238 239 ! Count the number of power_strain avec the previous check 240 ! or just remove the power_strain equal to zero 241 nstrain_tmp = 0 242 do idisp1=1,nstrain 243 if(power_strain_tmp(idisp1) > zero) then 244 nstrain_tmp = nstrain_tmp + 1 245 end if 246 end do 247 248 else 249 ndisp_tmp = ndisp 250 nstrain_tmp = nstrain 251 end if!end check 252 253 !init the values 254 polynomial_term%ndisp = ndisp_tmp 255 polynomial_term%nstrain = nstrain_tmp 256 polynomial_term%weight = weight 257 258 ABI_MALLOC(polynomial_term%atindx,(2,polynomial_term%ndisp)) 259 ABI_MALLOC(polynomial_term%direction,(polynomial_term%ndisp)) 260 ABI_MALLOC(polynomial_term%cell,(3,2,polynomial_term%ndisp)) 261 ABI_MALLOC(polynomial_term%power_disp,(polynomial_term%ndisp)) 262 ABI_MALLOC(polynomial_term%power_strain,(polynomial_term%nstrain)) 263 ABI_MALLOC(polynomial_term%strain,(polynomial_term%nstrain)) 264 265 !Transfert displacement 266 idisp2 = 0 267 do idisp1=1,ndisp 268 if(power_disp_tmp(idisp1) > zero)then 269 idisp2 = idisp2 + 1 270 polynomial_term%direction(idisp2) = direction(idisp1) 271 polynomial_term%power_disp(idisp2) = power_disp_tmp(idisp1) 272 polynomial_term%atindx(:,idisp2) = atindx(:,idisp1) 273 polynomial_term%cell(:,:,idisp2) = cell(:,:,idisp1) 274 polynomial_term%power_disp(idisp2) = power_disp_tmp(idisp1) 275 end if 276 end do 277 278 !Transfert strain 279 idisp2 = 0 280 do idisp1=1,nstrain 281 if(power_strain_tmp(idisp1) > zero)then 282 idisp2 = idisp2 + 1 283 polynomial_term%power_strain(idisp2) = power_strain_tmp(idisp1) 284 polynomial_term%strain(idisp2) = strain(idisp1) 285 end if 286 end do 287 288 end subroutine polynomial_term_init
m_polynomial_term/polynomial_term_type [ Types ]
[ Top ] [ m_polynomial_term ] [ Types ]
NAME
polynomial_term_type
FUNCTION
Datatype for a terms (displacements or strain) related to a polynomial coefficient
SOURCE
47 type, public :: polynomial_term_type 48 49 integer :: ndisp = 0 50 ! Number of displacement for this terms 51 ! 1 for (X_y-O_y)^3, 2 for (X_y-O_y)(X_x-O_y)^2... 52 53 integer :: nstrain = 0 54 ! Number of strain for this terms 55 56 integer,allocatable :: atindx(:,:) 57 ! atindx(2,ndisp) 58 ! Indexes of the atoms a and b in the unit cell 59 60 integer,allocatable :: cell(:,:,:) 61 ! cell(3,2,ndisp) 62 ! indexes of the cell of the atom a and b 63 64 integer,allocatable :: direction(:) 65 ! direction(ndisp) 66 ! direction of the displacement 67 68 integer,allocatable :: strain(:) 69 ! strain(nstrain) 70 ! strain 71 72 integer,allocatable :: power_disp(:) 73 ! power_disp(ndisp) 74 ! power of the displacement 2 (X_z-O_z)^2 or 1 for (X_y-O_y)^1 75 76 integer,allocatable :: power_strain(:) 77 ! power_strain(nstrain) 78 ! power of the strain 2 (\eta_1)^2 or 1 for (\eta_1)^1 79 80 real(dp) :: weight = zero 81 ! weight of the term 82 83 end type polynomial_term_type
m_polynomial_term/terms_compare [ Functions ]
[ Top ] [ m_polynomial_term ] [ Functions ]
NAME
equal
FUNCTION
Compare two polynomial_term_dot
INPUTS
t1<type(polynomial_term)> = datatype of the first term t2<type(polynomial_term)> = datatype of the second term
OUTPUT
res = logical
SOURCE
384 pure function terms_compare(t1,t2) result (res) 385 !Arguments ------------------------------------ 386 implicit none 387 388 !Arguments ------------------------------------ 389 type(polynomial_term_type), intent(in) :: t1,t2 390 logical :: res 391 !local 392 !variable 393 integer :: ia,idisp1,idisp2,mu 394 logical :: found 395 !array 396 integer :: blkval(2,t1%ndisp+t1%nstrain) 397 ! ************************************************************************* 398 res = .true. 399 blkval(:,:) = 0 400 if(t1%ndisp==t2%ndisp.and.t1%nstrain==t2%nstrain)then 401 ! Check strain 402 blkval(:,:) = 0 403 do idisp1=1,t1%nstrain 404 if(blkval(1,t1%ndisp+idisp1)==1)cycle!already found 405 do idisp2=1,t2%nstrain 406 if(blkval(2,t1%ndisp+idisp2)==1)cycle!already found 407 found = .false. 408 if(t1%strain(idisp1) == t2%strain(idisp2).and.& 409 & t1%power_strain(idisp1) == t2%power_strain(idisp2))then 410 found=.true. 411 end if 412 if(found)then 413 blkval(1,t1%ndisp+idisp1) = 1 414 blkval(2,t1%ndisp+idisp2) = 1 415 end if 416 end do 417 end do 418 if(any(blkval(:,t1%ndisp+1:t1%ndisp+t1%nstrain) == 0))then 419 res = .false. 420 return 421 end if 422 ! Check displacement 423 do idisp1=1,t1%ndisp 424 if(blkval(1,idisp1)==1)cycle!already found 425 do idisp2=1,t2%ndisp 426 if(blkval(2,idisp2)==1)cycle!already found 427 found = .false. 428 if(t1%atindx(1,idisp1) == t2%atindx(1,idisp2).and.& 429 & t1%atindx(2,idisp1) == t2%atindx(2,idisp2).and.& 430 & t1%direction(idisp1) == t2%direction(idisp2).and.& 431 & t1%power_disp(idisp1) == t2%power_disp(idisp2))then 432 found=.true. 433 do ia=1,2 434 do mu=1,3 435 if(t1%cell(mu,ia,idisp1) /= t2%cell(mu,ia,idisp2))then 436 found = .false. 437 cycle 438 end if 439 end do 440 end do 441 if(found)then 442 blkval(1,idisp1) = 1 443 blkval(2,idisp2) = 1 444 end if 445 end if 446 end do 447 end do 448 if(any(blkval(:,:)==0))res = .false. 449 else 450 res = .false. 451 end if 452 453 end function terms_compare