TABLE OF CONTENTS
- m_polynomial_coeff/check_irreducibility
- m_polynomial_coeff/coeffs_compare
- m_polynomial_coeff/coeffs_list_append
- m_polynomial_coeff/coeffs_list_conc
- m_polynomial_coeff/coeffs_list_conc_onsite
- m_polynomial_coeff/coeffs_list_copy
- m_polynomial_coeff/computeCombinationFromList
- m_polynomial_coeff/computeNorder
- m_polynomial_coeff/generateTermsFromList
- m_polynomial_coeff/getCoeffFromList
- m_polynomial_coeff/polynomial_coeff_broadcast
- m_polynomial_coeff/polynomial_coeff_evaluate
- m_polynomial_coeff/polynomial_coeff_free
- m_polynomial_coeff/polynomial_coeff_getEvenAnhaStrain
- m_polynomial_coeff/polynomial_coeff_getList
- m_polynomial_coeff/polynomial_coeff_getName
- m_polynomial_coeff/polynomial_coeff_getNorder
- m_polynomial_coeff/polynomial_coeff_getOrder1
- m_polynomial_coeff/polynomial_coeff_init
- m_polynomial_coeff/polynomial_coeff_list_free
- m_polynomial_coeff/polynomial_coeff_MPIrecv
- m_polynomial_coeff/polynomial_coeff_MPIsend
- m_polynomial_coeff/polynomial_coeff_setCoefficient
- m_polynomial_coeff/polynomial_coeff_setName
- m_polynomial_coeff/polynomial_coeff_type
- m_polynomial_coeff/reduce_zero_combinations
- m_polynomial_coeff/sort_combination
- m_polynomial_coeff/sort_combination_list
m_polynomial_coeff/check_irreducibility [ Functions ]
NAME
check_irreducibility
FUNCTION
checks irreducibility of combination of terms with respect to crystal symmetry
INPUTS
combination: combination of integers icoeff representing term A_x-B_x to check list_combination: list of combination against which combinatino will be checked list_symcoeff: list of coefficients containig information connecting the integers in combination and list_combination to bodies (A_x-B_x) of difference of atomic displacements list_symstr: list of strain with symmetry properties connecting icoeff>ncoeff_sym to strains of given direction ncoeff_sym: number of bodies (A_x-B_x) containing all symmetric equivalents nsym: number of symmetries in the crystal system ncombination: number of combinations -> defines size of list_combination ndisp: number of displacements that is maximum order of combinations index_irred: list of index of irreducible terms (that are non-zero) in list_combination
OUTPUT
logical :: irreducible -> TRUE: no other equal term exists in list_combination -> FALSE: a other equivalent term exists already SIDE_EFFECT: if combination is irreducible (that means irreducible = .FALSE.) the index of the irreduc!! irreducible coefficient is added to index_irred
SOURCE
4363 function check_irreducibility(combination,list_combination,list_symcoeff,list_symstr,ncoeff_sym,nsym,& 4364 & ncombination,ndisp,index_irred) result(irreducible) 4365 !Arguments ------------------------------------ 4366 implicit none 4367 4368 !Arguments ------------------------------------ 4369 !scalar 4370 integer,intent(in) :: ncoeff_sym,ncombination,ndisp,nsym 4371 integer,intent(inout),allocatable :: index_irred(:) 4372 logical :: irreducible !output 4373 !array 4374 integer,intent(inout) :: combination(ndisp),list_combination(ndisp,ncombination) 4375 integer,intent(in) :: list_symcoeff(6,ncoeff_sym,nsym) 4376 integer,intent(in) :: list_symstr(6,nsym,2) 4377 !local 4378 !variable 4379 integer :: icombi, isym 4380 !array 4381 integer :: combination_cmp_tmp(ndisp) 4382 integer,allocatable :: index_irred_tmp(:) 4383 ! ************************************************************************* 4384 ! sort input combinations 4385 !write(std_out,*) "DEBUG call sort_combination_lsit" 4386 call sort_combination_list(list_combination,size(list_combination,1),size(list_combination,2)) 4387 !write(std_out,*) "DEBUG sort_combination in check_irreducibility" 4388 call sort_combination(combination,size(combination)) 4389 4390 !Initialize variables 4391 icombi = 1 4392 irreducible = .TRUE. 4393 !do icombi=1,ncombination 4394 ! write(std_out,*) "DEBUG: list_combination i:",icombi,": ", list_combination(:,icombi) 4395 !enddo 4396 !MS DEBUG 4397 ! write(std_out,*) "DEBUG: shape(list_combination)", shape(list_combination(:,:)) 4398 !write(std_out,*) "DEBUG: combination: ", combination 4399 !write(std_out,*) "DEBUG: ncombination: ", ncombination 4400 4401 !Loop over non reducable combinations 4402 do while(icombi <= size(index_irred) .and. irreducible .eqv. .TRUE.) 4403 !if(all(list_combination(:,index_irred(icombi)) /= 0))then !.and. any(list_combination(:,i) <= ncoeff_symsym))then 4404 !If term j is equivalent to term i delete it 4405 if(all(list_combination(:,index_irred(icombi)) == combination(:)))then 4406 irreducible = .FALSE. 4407 return 4408 else !else loop over symmetries to find symmetry operation that projects term j on i 4409 isym = 2 4410 do while(isym <= nsym) 4411 !Get equivalent term indexes for symmetry isym 4412 call symcomb(combination, combination_cmp_tmp, isym) 4413 !Sort the new symmetric indexes for comparision 4414 call sort_combination(combination_cmp_tmp,size(combination_cmp_tmp)) 4415 !Compare. If equivalent break (term not irreducible 4416 if(all(list_combination(:,index_irred(icombi)) == combination_cmp_tmp(:)))then 4417 irreducible = .FALSE. 4418 return 4419 else 4420 isym = isym + 1 4421 endif 4422 enddo !isym 4423 endif ! all(list_combinaton...) 4424 !end if ! any(list_combination /= 0 ) 4425 icombi = icombi + 1 4426 end do !i=1,ncombination 4427 4428 !If no equivalent found put ncombination +1 to index of irreducible combis 4429 ABI_MALLOC(index_irred_tmp,(size(index_irred))) 4430 index_irred_tmp = index_irred 4431 ABI_FREE(index_irred) 4432 ABI_MALLOC(index_irred,(size(index_irred_tmp)+1)) 4433 index_irred(:size(index_irred_tmp)) = index_irred_tmp 4434 ABI_FREE(index_irred_tmp) 4435 index_irred(size(index_irred)) = ncombination + 1 4436 !write(std_out,*) "DEBUG: index_irred", index_irred 4437 return 4438 4439 contains 4440 subroutine symcomb(combination, combination_cmp_tmp, isym) 4441 integer, intent(inout) :: combination(:), combination_cmp_tmp(:), isym 4442 integer :: idisp, istrain 4443 do idisp=1,ndisp 4444 if(combination(idisp) /= 0 .and. combination(idisp) <= ncoeff_sym)then 4445 combination_cmp_tmp(idisp)=list_symcoeff(6,combination(idisp),isym) 4446 else if(combination(idisp) > ncoeff_sym)then 4447 istrain = combination(idisp) - ncoeff_sym 4448 combination_cmp_tmp(idisp)=list_symstr(istrain,isym,1) + ncoeff_sym 4449 else 4450 combination_cmp_tmp(idisp) = 0 4451 endif 4452 enddo 4453 end subroutine symcomb 4454 4455 end function check_irreducibility
m_polynomial_coeff/coeffs_compare [ Functions ]
NAME
equal
FUNCTION
INPUTS
OUTPUT
SOURCE
4026 function coeffs_compare(c1,c2) result (res) 4027 !Arguments ------------------------------------ 4028 implicit none 4029 4030 !Arguments ------------------------------------ 4031 type(polynomial_coeff_type), intent(in) :: c1,c2 4032 logical :: res 4033 !local 4034 !variable 4035 integer :: iterm1,iterm2 4036 !array 4037 !integer,allocatable :: blkval(:,:) 4038 integer :: blkval(2, max(c1%nterm,c2%nterm)) 4039 ! ************************************************************************* 4040 res = .false. 4041 blkval = 0 4042 do iterm1=1,c1%nterm 4043 if(blkval(1,iterm1)==1)cycle!already found 4044 do iterm2=1,c2%nterm 4045 if(blkval(2,iterm2)==1)cycle!already found 4046 if(c1%terms(iterm1)==c2%terms(iterm2)) then 4047 blkval(1,iterm1) = 1 4048 blkval(2,iterm2) = 1 4049 end if 4050 end do 4051 end do 4052 if(.not.any(blkval(:,:)==0))res = .true. 4053 4054 end function coeffs_compare
m_polynomial_coeff/coeffs_list_append [ Functions ]
NAME
coeff_list_append
FUNCTION
append on coeff to the list of coeffs: coeff_list
INPUTS
OUTPUT
SOURCE
4167 subroutine coeffs_list_append(coeff_list,coeff, check) 4168 !Arguments ------------------------------------ 4169 implicit none 4170 4171 !Arguments ------------------------------------ 4172 type(polynomial_coeff_type), allocatable, intent(inout) ::coeff_list(:) 4173 type(polynomial_coeff_type), intent(inout) :: coeff 4174 logical, intent(in) :: check 4175 !local variable 4176 type(polynomial_coeff_type) :: tmp(size(coeff_list)) 4177 integer :: n1, n2 4178 !array 4179 ! ************************************************************************* 4180 4181 n1=size(coeff_list) 4182 n2= n1+1 4183 4184 ! copy list1 to tmp 4185 tmp=coeff_list 4186 4187 ! allocate new list1 4188 call polynomial_coeff_list_free(coeff_list) 4189 ABI_MALLOC(coeff_list,(n2)) 4190 4191 call coeffs_list_copy(coeff_list, tmp) 4192 call polynomial_coeff_init(coeff%coefficient,coeff%nterm,coeff_list(n2),coeff%terms,& 4193 & coeff%name,check=check) 4194 end subroutine coeffs_list_append
m_polynomial_coeff/coeffs_list_conc [ Functions ]
NAME
coeff_list_conc
FUNCTION
Concatenate list1 and list2 of type polynomial_coeff and store it in list_out
INPUTS
OUTPUT
SOURCE
4072 function coeffs_list_conc(coeff_list1,coeff_list2) result (coeff_list_out) 4073 !Arguments ------------------------------------ 4074 implicit none 4075 4076 !Arguments ------------------------------------ 4077 type(polynomial_coeff_type), intent(in) :: coeff_list1(:),coeff_list2(:) 4078 type(polynomial_coeff_type) :: coeff_list_out(size(coeff_list1)+size(coeff_list2)) 4079 !local 4080 !variable 4081 integer :: ncoeff1,ncoeff2,ncoeff_out,i,j 4082 !array 4083 ! ************************************************************************* 4084 4085 !Get sizes of coeff_list1/2 4086 ncoeff1 = size(coeff_list1) 4087 ncoeff2 = size(coeff_list2) 4088 ncoeff_out = ncoeff1 + ncoeff2 4089 4090 ! ABI_MALLOC(coeff_list_out,(ncoeff_out)) 4091 do i=1,ncoeff_out 4092 if(i<=ncoeff1)then 4093 call polynomial_coeff_init(coeff_list1(i)%coefficient,coeff_list1(i)%nterm,coeff_list_out(i),coeff_list1(i)%terms,& 4094 & coeff_list1(i)%name,check=.TRUE.) 4095 else 4096 j=i-ncoeff1 4097 call polynomial_coeff_init(coeff_list2(j)%coefficient,coeff_list2(j)%nterm,coeff_list_out(i),coeff_list2(j)%terms,& 4098 & coeff_list2(j)%name,check=.TRUE.) 4099 endif 4100 enddo 4101 4102 end function coeffs_list_conc
m_polynomial_coeff/coeffs_list_conc_onsite [ Functions ]
NAME
coeff_list_conc_onsite
FUNCTION
Concatenate list1 and list2 of type polynomial_coeff and store it in list1
INPUTS
OUTPUT
SOURCE
4119 subroutine coeffs_list_conc_onsite(coeff_list1,coeff_list2) 4120 !Arguments ------------------------------------ 4121 implicit none 4122 4123 !Arguments ------------------------------------ 4124 type(polynomial_coeff_type), allocatable, intent(inout) :: coeff_list1(:) 4125 type(polynomial_coeff_type), intent(in) ::coeff_list2(:) 4126 !local 4127 !variable 4128 integer :: ncoeff1,ncoeff2,ncoeff_out 4129 !array 4130 type(polynomial_coeff_type), allocatable :: coeff_list_tmp(:) 4131 ! ************************************************************************* 4132 4133 ncoeff1=size(coeff_list1) 4134 ncoeff2=size(coeff_list2) 4135 ncoeff_out = ncoeff1+ncoeff2 4136 4137 ! copy list1 to tmp 4138 ABI_MALLOC(coeff_list_tmp,(ncoeff1)) 4139 coeff_list_tmp=coeff_list1 4140 4141 ! allocate new list1 4142 call polynomial_coeff_list_free(coeff_list1) 4143 ABI_MALLOC(coeff_list1,(ncoeff_out)) 4144 !coeff_list1=coeff_list_tmp + coeff_list2 4145 !coeff_list1 = coeffs_list_conc(coeff_list_tmp, coeff_list2) 4146 coeff_list1(1:ncoeff1) = coeff_list_tmp 4147 coeff_list1(ncoeff1+1:ncoeff_out) = coeff_list2 4148 call polynomial_coeff_list_free(coeff_list_tmp) 4149 4150 end subroutine coeffs_list_conc_onsite
m_polynomial_coeff/coeffs_list_copy [ Functions ]
NAME
coeff_list_copy
FUNCTION
Copy list1 to list2 of type polynomial_coeff
INPUTS
OUTPUT
SOURCE
4211 subroutine coeffs_list_copy(coeff_list_out,coeff_list_in) 4212 !Arguments ------------------------------------ 4213 implicit none 4214 4215 !Arguments ------------------------------------ 4216 type(polynomial_coeff_type), intent(in) :: coeff_list_in(:) 4217 type(polynomial_coeff_type),intent(out) :: coeff_list_out(:) 4218 !local 4219 !variable 4220 integer :: ncoeff_in,ncoeff_out,ii 4221 logical :: check 4222 character(len=500):: message 4223 !array 4224 ! ************************************************************************* 4225 4226 check = .true. 4227 4228 !Get size of coeff_lists 4229 ncoeff_in = size(coeff_list_in) 4230 ncoeff_out = size(coeff_list_out) 4231 4232 if(ncoeff_in > ncoeff_out)then 4233 write(message,'(a,a,a)')'The input list of polynomial_coefficients is larger',ch10,& 4234 & 'than the output list you want it assign to. Check size of lists.' 4235 ABI_ERROR(message) 4236 endif 4237 4238 !Copy input list into output lists 4239 do ii=1,ncoeff_in 4240 call polynomial_coeff_init(coeff_list_in(ii)%coefficient,coeff_list_in(ii)%nterm,& 4241 & coeff_list_out(ii),coeff_list_in(ii)%terms,coeff_list_in(ii)%name,check) 4242 enddo 4243 end subroutine coeffs_list_copy
m_polynomial_coeff/computeCombinationFromList [ Functions ]
NAME
computeCombinationFromList
FUNCTION
Recursive routine to compute the order N of a all the possible coefficient from the list list_symcoeff and list_symstr.
INPUTS
cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...) compatibleCoeffs(ncoeff+nstr,ncoeff+nstr) = array with the list of compatible coefficients 0 or 1 list_coeff(6,ncoeff_sym,nsym) = array with the list of the coefficients, for each coefficients (ncoeff_sym), we store the symmetrics(nsym) the 6th first dimensions are : 1 = direction of the IFC 2 = index of the atom number 1 (1=>natom) 3 = index of the atom number 2 (1=>natom) 4 = indexes of the cell of the second atom (the atom number 1 is always in the cell 0 0 0) 5 = weight of the term (-1 or 1) 6 = indexes of the symmetric list_str(nstr_sym,nsym) = array with the list of the strain and the symmetrics index_coeff_in(power_disp-1) = list of previous coefficients computed (start with 0) icoeff = current indexes of the combination (start with 1) max_power_strain = maximum order of the strain of the strain phonon coupling nmodel_tot = current number of combination already computed (start with 0) natom = number of atoms in the unit cell ncoeff = number of coefficient for related to the atomic displacment into list_symcoeff nstr = number of coefficient for related to the strain into list_symstr nmodel = number of maximum models nrpt = number of cell nsym = number of symmetries in the system For example, the sum of all the term like (Sr_y-O_y)^odd, are 0 by symetrie in cubic system. Here, we build a list with: 0 this term is not allowed for odd 1 this term is allowed for odd power_disp = initial power_disp to be computed (can be < power_disp_min, this routine will skip the first power_disp) power_disp_min = minimal power_disp to be computed power_disp_max = maximum power_disp to be computed symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...) nbody = optional, number of body for the coefficients, for example: 0 => all the terms 1 => only (Sr_x-T_y)^power_disp and (Sr_x-T_y)^power_disp\eta^power_disp ... compute = logical, optional: TRUE if we store the coefficients FALSE just to count the number of coefficient anharmstr = logical, optional : TRUE, the anharmonic strain are computed FALSE, (default) the anharmonic strain are not computed distributed = logical, optional : True, the coefficients will be distributed on the CPU only_odd_power = logical, optional : if TRUE return only odd power only_even_power= logical, optional : if TRUe return only even power
OUTPUT
icoeff = current indexes of the cofficients (start we 1) nmodel_tot = current number of coefficients already computed (start we 0) list_combination = list of the possible combination of coefficients
SOURCE
3131 recursive subroutine computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,& 3132 & index_coeff_in,list_combination,icoeff,max_power_strain,& 3133 & natom,ncoeff,ncoeff_sym,iirred_comb,nirred_comb,nstr,nmodel,nrpt,nsym,power_disp,power_disp_min,& 3134 & power_disp_max,symbols,comm,nbody,only_odd_power,only_even_power,& 3135 & compute,anharmstr,spcoupling,disp) 3136 3137 implicit none 3138 3139 !Arguments --------------------------------------------- 3140 !scalar 3141 integer,intent(in) :: natom,ncoeff,ncoeff_sym,power_disp,power_disp_min,power_disp_max 3142 integer,intent(in) :: max_power_strain,nmodel,nsym,nrpt,nstr,comm,icoeff 3143 integer,intent(inout) :: nirred_comb,iirred_comb 3144 logical,optional,intent(in) :: compute,anharmstr,spcoupling,disp 3145 integer,optional,intent(in) :: nbody 3146 logical,optional,intent(in) :: only_odd_power,only_even_power 3147 !arrays 3148 integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff_sym+nstr,ncoeff_sym+nstr) 3149 integer,intent(in) :: list_coeff(6,ncoeff_sym,nsym),list_str(nstr,nsym,2) 3150 integer,intent(in) :: index_coeff_in(power_disp-1) 3151 integer,intent(out) :: list_combination(power_disp_max,nirred_comb) 3152 character(len=5),intent(in) :: symbols(natom) 3153 !Local variables --------------------------------------- 3154 !scalar 3155 integer :: icoeff1,icoeff2,nbody_in,ii,jj,nbody_count 3156 integer :: ndisp_out,nstrain 3157 logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling 3158 logical :: need_only_odd_power,need_only_even_power,compute_sym,need_disp 3159 !arrays 3160 integer :: powers(power_disp) 3161 integer,allocatable :: index_coeff(:) 3162 ! ************************************************************************* 3163 3164 !Set the inputs 3165 need_compute = .TRUE. 3166 need_anharmstr = .TRUE. 3167 need_spcoupling = .TRUE. 3168 need_disp = .TRUE. 3169 need_only_odd_power = .FALSE. 3170 need_only_even_power = .FALSE. 3171 compute_sym = .FALSE. !Never compute the symmetric combinations here 3172 nbody_in = 0 !all kind of terms 3173 if(present(compute)) need_compute = compute 3174 if(present(nbody)) nbody_in = nbody 3175 if(present(anharmstr)) need_anharmstr = anharmstr 3176 if(present(spcoupling)) need_spcoupling = spcoupling 3177 if(present(disp)) need_disp = disp 3178 if(present(only_odd_power)) need_only_odd_power = only_odd_power 3179 if(present(only_even_power)) need_only_even_power = only_even_power 3180 3181 if(power_disp <= power_disp_max)then 3182 3183 ! Initialisation of variables 3184 ABI_MALLOC(index_coeff,(power_disp)) 3185 index_coeff(1:power_disp-1) = index_coeff_in(:) 3186 ! write(std_out,*) "DEBUG index_coff in recursive CCL: ", index_coeff 3187 ! Loop over ncoeff+nstr 3188 do icoeff1=icoeff,ncoeff+nstr 3189 3190 ! Reset the flag compatible and possible 3191 compatible = .TRUE. 3192 possible = .TRUE. 3193 3194 ! If the power_disp is one, we need to set icoeff to icoeff1 3195 if(power_disp==1) then 3196 if(icoeff1<=ncoeff .and. compatibleCoeffs(icoeff,icoeff1)==0)then 3197 ! is_displacement and compatible 3198 compatible = .FALSE. 3199 end if 3200 end if 3201 ! If the distance between the 2 coefficients is superior than the cut-off, we cycle. 3202 do icoeff2=1,power_disp-1 3203 ! write(std_out,*) "icoeff1: ", icoeff1 3204 ! write(std_out,*) "icoeff2: ", icoeff2, "index_icoeff2: ", index_coeff(icoeff2) 3205 if(icoeff1 <= ncoeff .and. index_coeff(icoeff2) <=ncoeff)then 3206 if(compatibleCoeffs(index_coeff(icoeff2),icoeff1)==0)then 3207 compatible = .FALSE. 3208 end if 3209 endif 3210 end do 3211 3212 if (.not.compatible) cycle !The distance is not compatible 3213 3214 ! Set the index of the new coeff in the list 3215 index_coeff(power_disp) = icoeff1 3216 ! Do some checks 3217 ! ------------- 3218 ! 1-Check if the coefficient is full anharmonic strain and if we need to compute it 3219 if(all(index_coeff > ncoeff))then 3220 compatible = (need_anharmstr .or. need_spcoupling) 3221 possible = need_anharmstr 3222 end if 3223 ! 2-Check if the coefficient is strain-coupling and if we need to compute it 3224 if(any(index_coeff <= ncoeff) .and. any(index_coeff > ncoeff))then 3225 ! write(std_out,*) "index_coeff", index_coeff,"need_spcoupling",need_spcoupling 3226 possible = need_spcoupling 3227 compatible = need_spcoupling 3228 if(count(index_coeff > ncoeff) > max_power_strain)then 3229 possible = .false. 3230 compatible = .false. 3231 end if 3232 end if 3233 ! 3-Check if the coefficient is only disp and if we need to compute it 3234 if(all(index_coeff <= ncoeff))then 3235 ! write(std_out,*) "index_coeff", index_coeff,"need_dis",need_disp 3236 compatible = (need_disp .or. need_spcoupling) 3237 possible = need_disp 3238 end if 3239 ! 4-Count number of Strain and number of displacements for compute symmetric terms 3240 nstrain = 0 3241 ndisp_out = 0 3242 do ii=1,power_disp 3243 if(index_coeff(ii) > 0 .and. index_coeff(ii) <= ncoeff)then 3244 ndisp_out = ndisp_out + 1 3245 else 3246 nstrain = nstrain +1 3247 index_coeff(ii) = index_coeff(ii) - ncoeff + ncoeff_sym 3248 end if 3249 end do 3250 3251 if(power_disp >= power_disp_min) then 3252 3253 ! count the number of body 3254 powers(:) = 1 3255 do ii=1,power_disp 3256 do jj=ii+1,power_disp 3257 if (powers(jj) == 0) cycle 3258 if(index_coeff(ii)==index_coeff(jj))then 3259 powers(ii) = powers(ii) + 1 3260 powers(jj) = 0 3261 end if 3262 end do 3263 end do 3264 nbody_count = count(powers /= 0) 3265 ! write(std_out,*) "powers: ", powers 3266 ! write(std_out,*) "nbody_count: ", nbody_count 3267 ! check the only_odd and only_even flags 3268 if(any(mod(powers(1:power_disp),2) /=0) .and. need_only_even_power) then 3269 possible = .false. 3270 end if 3271 if(any(mod(powers(1:power_disp),2) ==0) .and. need_only_odd_power)then 3272 possible = .false. 3273 end if 3274 ! Check the nbody flag 3275 if(nbody_in /= 0)then 3276 if(power_disp-count(powers==0) > nbody_in) then 3277 possible = .false. 3278 compatible = .false. 3279 end if 3280 end if 3281 3282 if(possible) then 3283 ! increase coefficients and set it 3284 ! nmodel_tot = nmodel_tot + 1 3285 ! if(need_compute)then 3286 ! list_combination(1:power_disp,nmodel_tot) = index_coeff 3287 ! end if 3288 !nmodel_tot_test = 0 3289 ! !Start from second symmetry in Symmetric Combinations 3290 ! isym_in_test = 2 3291 ! idisp_in_test = power_disp 3292 ! ndisp_test = power_disp 3293 ! index_coeff_tmp = index_coeff 3294 !Count anharmonic strain terms 3295 if(ndisp_out == 0 .and. nstrain > 0)then 3296 nirred_comb = nirred_comb +1 3297 iirred_comb = iirred_comb +1 3298 if(need_compute)then 3299 list_combination(1:power_disp,iirred_comb) = index_coeff 3300 endif 3301 else !Else counst symmetric terms of atomic displacement (pure disp or disp/strain) 3302 !Store index for each combination of irreducible terms to later parallely compute symmetric combinations 3303 nirred_comb = nirred_comb +1 3304 iirred_comb = iirred_comb +1 3305 !write(std_out,*) "DEBUG index_coeff: ", index_coeff 3306 !write(std_out,*) "DEBUG ndisp_out: ", ndisp_out 3307 !write(std_out,*) "DEBUG nstrain: ", nstrain 3308 !write(std_out,*) "DEBUG nmodel_start: ", nmodel_start 3309 if(need_compute)then 3310 list_combination(:,iirred_comb) = 0 3311 list_combination(:ndisp_out+nstrain,iirred_comb) = index_coeff 3312 endif 3313 end if !ndisp_out == 0 .and.n nstrain >0 3314 end if!possible 3315 end if!end if power_disp < power_disp_min 3316 3317 !Change back to irreducible terms ncoeff_limit 3318 do ii=1,power_disp 3319 if(index_coeff(ii) > ncoeff_sym)then 3320 index_coeff(ii) = index_coeff(ii) + ncoeff - ncoeff_sym 3321 end if 3322 end do 3323 3324 ! If the model is still compatbile with the input flags, we continue. 3325 if(compatible)then 3326 call computeCombinationFromList(cell,compatibleCoeffs,list_coeff,list_str,& 3327 & index_coeff,list_combination,icoeff1,max_power_strain,& 3328 & natom,ncoeff,ncoeff_sym,iirred_comb,nirred_comb,nstr,nmodel,nrpt,nsym,power_disp+1,& 3329 & power_disp_min,power_disp_max,symbols,comm,nbody=nbody_in,& 3330 & compute=need_compute,anharmstr=need_anharmstr,& 3331 & spcoupling=need_spcoupling,only_odd_power=need_only_odd_power,& 3332 & only_even_power=need_only_even_power,disp=need_disp) 3333 end if 3334 end do 3335 ABI_FREE(index_coeff) 3336 end if 3337 3338 end subroutine computeCombinationFromList
m_polynomial_coeff/computeNorder [ Functions ]
NAME
computeNorder
FUNCTION
Recursive routine to compute the order N of a all the possible coefficient from the list list_symcoeff and list_symstr.
INPUTS
cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...) compatibleCoeffs(ncoeff+nstr,ncoeff+nstr) = array with the list of compatible coefficients 0 or 1 list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients, for each coefficients (ncoeff_sym), we store the symmetrics(nsym) the 6th first dimensions are : 1 = direction of the IFC 2 = index of the atom number 1 (1=>natom) 3 = index of the atom number 2 (1=>natom) 4 = indexes of the cell of the second atom (the atom number 1 is always in the cell 0 0 0) 5 = weight of the term (-1 or 1) 6 = indexes of the symmetric list_symstr(nstr_sym,nsym) = array with the list of the strain and the symmetrics index_coeff_in(power_disp-1) = list of previous coefficients computed (start with 0) icoeff = current indexes of the cofficients (start we 1) icoeff_tot = current number of coefficients already computed (start we 0) natom = number of atoms in the unit cell nstr = number of coefficient for related to the atomic displacment into list_symcoeff nstr = number of coefficient for related to the strain into list_symstr ncoeff_out = number of maximum coefficients nrpt = number of cell nsym = number of symmetries in the system power_disp = initial power_disp to be computed (can be < power_disp_min, this routine will skip the firts power_disp) power_disp_min = minimal power_disp to be computed power_disp_max = maximum power_disp to be computed symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...) nbody = optional, number of body for the coefficients, for example: 0 => all the terms 1 => only (Sr_x-T_y)^power_disp and (Sr_x-T_y)^power_disp\eta^power_disp ... compute = logical, optional: TRUE if we store the coefficients FALSE just to count the number of coefficient anharmstr = logical, optional : TRUE, the anharmonic strain are computed FALSE, (default) the anharmonic strain are not computed distributed = logical, optional : True, the coefficients will be distributed on the CPU
OUTPUT
icoeff = current indexes of the cofficients (start we 1) icoeff_tot = current number of coefficients already computed (start we 0) polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with the polynomial_coeff
SOURCE
2883 recursive subroutine computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,& 2884 & index_coeff_in,icoeff,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,& 2885 & nrpt,nsym,power_disp,power_disp_min,power_disp_max,symbols,nbody,& 2886 & compute,anharmstr,spcoupling,distributed) 2887 2888 implicit none 2889 2890 !Arguments --------------------------------------------- 2891 !scalar 2892 integer,intent(in) :: natom,ncoeff,power_disp,power_disp_min,power_disp_max,ncoeff_out,nsym,nrpt,nstr,icoeff 2893 integer,intent(inout) :: icoeff_tot 2894 logical,optional,intent(in) :: compute,anharmstr,spcoupling,distributed 2895 integer,optional,intent(in) :: nbody 2896 !arrays 2897 integer,intent(in) :: cell(3,nrpt),compatibleCoeffs(ncoeff+nstr,ncoeff+nstr) 2898 integer,intent(in) :: list_coeff(6,ncoeff,nsym),list_str(nstr,nsym,2) 2899 integer,intent(in) :: index_coeff_in(power_disp-1) 2900 type(polynomial_coeff_type),intent(inout) :: coeffs_out(ncoeff_out) 2901 character(len=5),intent(in) :: symbols(natom) 2902 !Local variables --------------------------------------- 2903 !scalar 2904 integer :: ia,ib,ii,icoeff1,icoeff_tmp 2905 integer :: iterm,nbody_in,ncoeff_max,pa,pb 2906 integer :: ndisp_max,nterm_max 2907 real(dp):: coefficient 2908 logical :: need_compute,compatible,possible,need_anharmstr,need_spcoupling,need_distributed 2909 !arrays 2910 integer,allocatable :: index_coeff(:) 2911 character(len=200):: name 2912 type(polynomial_term_type),dimension(:),allocatable :: terms 2913 type(polynomial_coeff_type),allocatable :: coeffs_tmp(:) 2914 ! ************************************************************************* 2915 2916 !Set the inputs 2917 need_compute = .TRUE. 2918 need_anharmstr = .TRUE. 2919 need_spcoupling = .TRUE. 2920 need_distributed = .FALSE. 2921 nbody_in = 0 !all kind of terms 2922 if(present(compute)) need_compute = compute 2923 if(present(nbody)) nbody_in = nbody 2924 if(present(anharmstr)) need_anharmstr = anharmstr 2925 if(present(spcoupling)) need_spcoupling = spcoupling 2926 if(present(distributed)) need_distributed = distributed 2927 if(power_disp <= power_disp_max)then 2928 2929 ! Initialisation of variables 2930 nterm_max = nsym 2931 ncoeff_max = (ncoeff+nstr) 2932 ndisp_max = power_disp 2933 icoeff_tmp = 0 2934 ABI_MALLOC(coeffs_tmp,(ncoeff_max)) 2935 ABI_MALLOC(terms,(nterm_max)) 2936 ABI_MALLOC(index_coeff,(power_disp)) 2937 2938 index_coeff(1:power_disp-1) = index_coeff_in(:) 2939 2940 do icoeff1=icoeff,ncoeff+nstr 2941 ! If the distance between the 2 coefficients is superior than the cut-off, 2942 ! we cycle 2943 ! If the power_disp is one check if icoeff1 is compatible with itself 2944 if(power_disp==1) then 2945 if(icoeff1 <= ncoeff .and. compatibleCoeffs(icoeff1,icoeff1)==0)then 2946 cycle 2947 end if 2948 end if 2949 if(compatibleCoeffs(icoeff,icoeff1)==0) cycle 2950 2951 ! Reset the flag compatible and possible 2952 compatible = .TRUE. 2953 possible = .TRUE. 2954 2955 index_coeff(power_disp) = icoeff1 2956 iterm = 0 2957 coefficient = one 2958 2959 if(power_disp >= power_disp_min) then 2960 block 2961 logical :: reverse(ndisp_max) 2962 reverse(:) = .False. 2963 call generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,& 2964 & ndisp_max,nrpt,nstr,nsym,iterm,terms, reverse=reverse) 2965 end block 2966 2967 if(iterm > 0)then 2968 ! Do some checks 2969 ! ------------- 2970 ! 1-Check if the coefficient is full anharmonic strain and if we need to compute it 2971 if(terms(1)%ndisp == 0)then 2972 compatible = (need_anharmstr .or. need_spcoupling) 2973 possible = need_anharmstr 2974 end if 2975 ! 1-Check if the coefficient is strain-coupling and if we need to compute it 2976 if(terms(1)%nstrain > 0.and.terms(1)%ndisp > 0)then 2977 possible = need_spcoupling 2978 compatible = need_spcoupling 2979 end if 2980 2981 ! ------------ 2982 ! 2-Check if this terms is compatible with nbody 2983 if(nbody_in > 0)then 2984 pa = 1 ; pb = 1 2985 ia = 0 ; ib = 0 2986 ! Count the number of terms and the power_disp 2987 do ii=1,terms(1)%ndisp 2988 if(terms(1)%nstrain > 0) then 2989 pb = pb*terms(1)%power_disp(ii) 2990 ib = ib + 1 2991 else 2992 pa = pa*terms(1)%power_disp(ii) 2993 ia = ia + 1 2994 end if 2995 end do 2996 if(ia <= nbody_in)then 2997 if(ia==nbody_in.and.abs(mod(pa,2)) < tol16)then 2998 if(ib==0)then 2999 compatible = .FALSE. 3000 possible = .TRUE. 3001 else if (ib==nbody_in.and.abs(mod(pb,2)) < tol16) then 3002 compatible = .FALSE. 3003 possible = .TRUE. 3004 else 3005 possible = .FALSE. 3006 compatible = .FALSE. 3007 end if 3008 else 3009 possible = .FALSE. 3010 compatible = .FALSE. 3011 end if 3012 else 3013 compatible = .FALSE. 3014 possible = .FALSE. 3015 end if 3016 end if 3017 3018 if(possible)then 3019 ! increase coefficients and set it 3020 icoeff_tmp = icoeff_tmp + 1 3021 icoeff_tot = icoeff_tot + 1 3022 call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp),& 3023 & terms(1:iterm),check=.true.) 3024 3025 end if 3026 end if 3027 3028 ! Deallocate the terms 3029 do iterm=1,nterm_max 3030 call polynomial_term_free(terms(iterm)) 3031 end do 3032 3033 end if!end if power_disp < power_disp_min 3034 3035 if(compatible)then 3036 call computeNorder(cell,coeffs_out,compatibleCoeffs,list_coeff,list_str,index_coeff,& 3037 & icoeff1,icoeff_tot,natom,ncoeff,nstr,ncoeff_out,nrpt,nsym,power_disp+1,& 3038 & power_disp_min,power_disp_max,symbols,nbody=nbody_in,compute=need_compute,& 3039 & anharmstr=need_anharmstr,spcoupling=need_spcoupling) 3040 end if 3041 end do 3042 3043 ABI_FREE(terms) 3044 ABI_FREE(index_coeff) 3045 3046 ! Transfer in the final array 3047 icoeff1 = 0 3048 do icoeff_tmp=1,ncoeff_max 3049 if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then 3050 ! Increase icoeff and fill the coeffs_out array 3051 icoeff_tot = icoeff_tot + 1 3052 if(need_compute)then 3053 name = '' 3054 ! Get the name of this coefficient 3055 call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.) 3056 call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,& 3057 & coeffs_out(icoeff_tot),coeffs_tmp(icoeff_tmp)%terms,& 3058 & name=name) 3059 end if 3060 end if 3061 end do 3062 ! Deallocation 3063 do icoeff1=1,ncoeff_max 3064 call polynomial_coeff_free(coeffs_tmp(icoeff1)) 3065 end do 3066 ABI_FREE(coeffs_tmp) 3067 end if 3068 3069 end subroutine computeNorder
m_polynomial_coeff/generateTermsFromList [ Functions ]
NAME
generateTermsFromList
FUNCTION
Compute for a given list of index the correspondig set of terms
INPUTS
cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...) index_coeff_in(ndisp) = list of coefficients to be computed list_symcoeff(6,ncoeff,nsym) = array with the list of the coefficients, for each coefficients (ncoeff_sym), we store the symmetrics(nsym) the 6th first dimensions are : 1 = direction of the IFC 2 = index of the atom number 1 (1=>natom) 3 = index of the atom number 2 (1=>natom) 4 = indexes of the cell of the second atom (the atom number 1 is always in the cell 0 0 0) 5 = weight of the term (-1 or 1) 6 = indexes of the symmetric list_symstr(nstr,nsym) = array with the list of the strain and the symmetrics ncoeff = number of maximum coefficients in the list_symcoeff ndisp = number of maximum diplacement (phonon + strain) nrpt = number of cell nsym = number of symmetries in the system
OUTPUT
terms<(type(polynomial_term_type)>(nterm) = list of terms nterm = number of ouput terms
SOURCE
3645 subroutine generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,ndisp_max,& 3646 & nrpt,nstr,nsym,nterm,terms, reverse) 3647 3648 implicit none 3649 3650 !Arguments ------------------------------------ 3651 !scalar 3652 integer,intent(in) :: ndisp_max,ncoeff,nrpt,nstr,nsym 3653 integer,intent(out):: nterm 3654 !arrays 3655 integer,intent(in) :: index_coeff(ndisp_max) 3656 integer,intent(in) :: cell(3,nrpt),list_coeff(6,ncoeff,nsym) 3657 integer,intent(in) :: list_str(nstr,nsym,2) 3658 logical, intent(in) :: reverse(ndisp_max) 3659 type(polynomial_term_type),intent(out) :: terms(nsym) 3660 !Local variables------------------------------- 3661 !scalar 3662 integer :: ia,ib,icoeff_str,idisp,irpt 3663 integer :: isym,ndisp,nstrain,mu 3664 real(dp):: weight 3665 !arrays 3666 integer :: atindx(2,ndisp_max),cells(3,2,ndisp_max),dir_int(ndisp_max),strain(ndisp_max) 3667 integer :: power_disps(ndisp_max),power_strain(ndisp_max) 3668 3669 ! ************************************************************************* 3670 nterm = 0 3671 !Loop over symetries 3672 do isym=1,nsym 3673 !Treat this coeff 3674 weight = 1 3675 ndisp = 0 3676 nstrain = 0 3677 do idisp=1,ndisp_max 3678 ! Get index of this displacement term 3679 ! Check if the index is not zero 3680 if(index_coeff(idisp)==0)cycle 3681 if(index_coeff(idisp)<=ncoeff)then 3682 ndisp = ndisp + 1 3683 mu = list_coeff(1,index_coeff(idisp),isym) 3684 if( reverse(idisp)) then 3685 ia = list_coeff(3,index_coeff(idisp),isym) 3686 ib = list_coeff(2,index_coeff(idisp),isym) 3687 irpt = list_coeff(4,index_coeff(idisp),isym) 3688 irpt = find_opposite_irpt(cell, irpt) 3689 weight = -weight*list_coeff(5,index_coeff(idisp),isym) 3690 else 3691 ia = list_coeff(2,index_coeff(idisp),isym) 3692 ib = list_coeff(3,index_coeff(idisp),isym) 3693 irpt = list_coeff(4,index_coeff(idisp),isym) 3694 weight = weight*list_coeff(5,index_coeff(idisp),isym) 3695 end if 3696 ! Fill First term arrays 3697 atindx(1,idisp) = ia; atindx(2,idisp) = ib; 3698 dir_int(idisp) = mu 3699 power_disps(idisp) = 1 3700 cells(:,1,idisp) = (/0,0,0/) 3701 cells(:,2,idisp) = cell(:,irpt) 3702 else 3703 nstrain = nstrain + 1 3704 icoeff_str = index_coeff(idisp)-ncoeff 3705 strain(nstrain) = list_str(icoeff_str,isym,1) 3706 power_strain(nstrain) = 1 3707 weight = weight*list_str(icoeff_str,isym,2) 3708 end if 3709 end do 3710 nterm = nterm + 1 3711 call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(nterm),power_disps,& 3712 & power_strain,strain,weight,check=.true.) 3713 end do!end do sym 3714 3715 3716 end subroutine generateTermsFromList
m_polynomial_coeff/getCoeffFromList [ Functions ]
NAME
getCoeffFromList
FUNCTION
get the index of a coefficient into the list_coeff
INPUTS
list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients, for each coefficients (ncoeff_sym), we store the symmetrics(nsym) the 6th first dimensions are : 1 = direction of the IFC 2 = index of the atom number 1 (1=>natom) 3 = index of the atom number 2 (1=>natom) 4 = indexes of the cell of the second atom (the atom number 1 is always in the cell 0 0 0) 5 = weight of the term (-1 or 1) 6 = indexes of the symmetric ia = index of the atom 1 ib = index of the atom 1 irpt = indexes of the cell of the second atom mu = direction of the IFC ncoeff = number of total coefficients in the list
OUTPUT
coeff = index of the coefficient
SOURCE
3579 function getCoeffFromList(list_coeff,ia,ib,irpt,mu,ncoeff) result(coeff) 3580 3581 implicit none 3582 3583 !Arguments ------------------------------------ 3584 !scalar 3585 integer,intent(in) :: ia,ib,irpt,mu,ncoeff 3586 integer :: coeff 3587 !arrays 3588 integer,intent(in) :: list_coeff(6,ncoeff) 3589 !Local variables------------------------------- 3590 !scalar 3591 integer :: icoeff 3592 !arrays 3593 3594 ! ************************************************************************* 3595 coeff = 0 3596 !write(std_out,*) "DEBUG shape inside FromList: ", shape(list_coeff) 3597 do icoeff = 1,ncoeff 3598 if(mu==list_coeff(1,icoeff).and.& 3599 & ia==list_coeff(2,icoeff).and.& 3600 & ib==list_coeff(3,icoeff).and.& 3601 & irpt==list_coeff(4,icoeff))then!.and.& 3602 !& abs(weight-list_coeff(5,icoeff)) < tol16) then 3603 coeff = icoeff 3604 exit 3605 end if 3606 end do 3607 3608 end function getCoeffFromList
m_polynomial_coeff/polynomial_coeff_broadcast [ Functions ]
NAME
polynomial_coeff_broadcast
FUNCTION
MPI broadcast polynomial_coefficent datatype
INPUTS
source = rank of source comm = MPI communicator
SIDE EFFECTS
coefficients<type(polynomial_coefficent_type)>= Input if node is source, other nodes returns with a completely initialized instance.
SOURCE
606 subroutine polynomial_coeff_broadcast(coefficients, source, comm) 607 608 implicit none 609 610 !Arguments ------------------------------------ 611 !array 612 type(polynomial_coeff_type),intent(inout) :: coefficients 613 integer, intent(in) :: source,comm 614 615 !Local variables------------------------------- 616 !scalars 617 integer :: ierr,ii 618 !arrays 619 620 ! ************************************************************************* 621 622 623 if (xmpi_comm_size(comm) == 1) return 624 625 ! Free the output 626 if (xmpi_comm_rank(comm) /= source) then 627 call polynomial_coeff_free(coefficients) 628 end if 629 630 ! Transmit variables 631 call xmpi_bcast(coefficients%name, source, comm, ierr) 632 call xmpi_bcast(coefficients%nterm, source, comm, ierr) 633 call xmpi_bcast(coefficients%coefficient, source, comm, ierr) 634 635 !Allocate arrays on the other nodes. 636 if (xmpi_comm_rank(comm) /= source) then 637 ABI_MALLOC(coefficients%terms,(coefficients%nterm)) 638 do ii=1,coefficients%nterm 639 call polynomial_term_free(coefficients%terms(ii)) 640 end do 641 end if 642 ! Set the number of term on each node (needed for allocations of array) 643 do ii = 1,coefficients%nterm 644 call xmpi_bcast(coefficients%terms(ii)%ndisp, source, comm, ierr) 645 call xmpi_bcast(coefficients%terms(ii)%nstrain, source, comm, ierr) 646 end do 647 648 ! Allocate arrays on the other nodes 649 if (xmpi_comm_rank(comm) /= source) then 650 do ii = 1,coefficients%nterm 651 ABI_MALLOC(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp)) 652 coefficients%terms(ii)%atindx = 0 653 ABI_MALLOC(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp)) 654 ABI_MALLOC(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp)) 655 ABI_MALLOC(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp)) 656 ABI_MALLOC(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain)) 657 ABI_MALLOC(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain)) 658 end do 659 end if 660 661 ! Transfert value 662 do ii = 1,coefficients%nterm 663 call xmpi_bcast(coefficients%terms(ii)%weight, source, comm, ierr) 664 call xmpi_bcast(coefficients%terms(ii)%atindx, source, comm, ierr) 665 call xmpi_bcast(coefficients%terms(ii)%direction, source, comm, ierr) 666 call xmpi_bcast(coefficients%terms(ii)%cell, source, comm, ierr) 667 call xmpi_bcast(coefficients%terms(ii)%power_disp, source, comm, ierr) 668 call xmpi_bcast(coefficients%terms(ii)%power_strain, source, comm, ierr) 669 call xmpi_bcast(coefficients%terms(ii)%strain, source, comm, ierr) 670 end do 671 672 673 end subroutine polynomial_coeff_broadcast
m_polynomial_coeff/polynomial_coeff_evaluate [ Functions ]
NAME
polynomial_coeff_evaluate
FUNCTION
Compute the energy related to the coefficients from fitted polynome
INPUTS
coefficients(ncoeff)<type(polynomial_coeff_type)> = list of coefficients disp(3,natom_sc) = atomics displacement between configuration and the reference natom_sc = number of atoms in the supercell natom_uc = number of atoms in the unit cell ncoeff = number of coefficients sc_size(3) = size of the supercell (2 2 2 for example) strain(6) = strain between configuration and the reference cells(ncell) = number of the cells into the supercell (1,2,3,4,5) ncell = total number of cell to treat by this cpu index_cells(3,ncell) = indexes of the cells into supercell (-1 -1 -1 ,...,1 1 1) comm=MPI communicator
OUTPUT
energy = contribution to the energy energy_coeff(ncoeff) = energy contribution of each anharmonic term fcart(3,natom) = contribution to the forces strten(6) = contribution to the stress tensor
SOURCE
1018 subroutine polynomial_coeff_evaluate(coefficients,disp,energy,energy_coeff,fcart,natom_sc,natom_uc,ncoeff,sc_size,& 1019 & strain,strten,ncell,index_cells,comm,filename) 1020 1021 !Arguments ------------------------------------ 1022 ! scalar 1023 real(dp),intent(out):: energy 1024 integer, intent(in) :: ncell,ncoeff,natom_sc,natom_uc 1025 integer, intent(in) :: comm 1026 character(len=fnlen),optional,intent(in) :: filename 1027 ! array 1028 real(dp),intent(out):: strten(6) 1029 real(dp),intent(in) :: strain(6) 1030 real(dp),intent(out):: fcart(3,natom_sc) 1031 real(dp),intent(in) :: disp(3,natom_sc) 1032 real(dp),optional,intent(out):: energy_coeff(ncoeff) 1033 integer,intent(in) :: index_cells(4,ncell) 1034 integer,intent(in) :: sc_size(3) 1035 type(polynomial_coeff_type),intent(in) :: coefficients(ncoeff) 1036 !Local variables------------------------------- 1037 ! scalar 1038 integer :: i1,i2,i3,ia1,ib1,ia2,ib2,idir1,idir2,ierr,ii 1039 integer :: icoeff,iterm,idisp1,idisp2,idisp1_strain,idisp2_strain,icell,ndisp 1040 integer :: nstrain,ndisp_tot,power_disp,power_strain,unit_out 1041 real(dp):: coeff,disp1,disp2,tmp1,tmp2,tmp3,weight 1042 logical :: file_opened 1043 ! array 1044 integer :: cell_atoma1(3),cell_atoma2(3) 1045 integer :: cell_atomb1(3),cell_atomb2(3) 1046 character(len=500) :: msg 1047 character(len=fnlen) :: name_file 1048 ! ************************************************************************* 1049 1050 ! Check 1051 if (any(sc_size <= 0)) then 1052 write(msg,'(a,a)')' No supercell found for getEnergy' 1053 ABI_ERROR(msg) 1054 end if 1055 1056 if(present(filename)) name_file = filename 1057 1058 ! Initialisation of variables 1059 energy = zero 1060 fcart(:,:) = zero 1061 strten(:) = zero 1062 energy_coeff(:) = zero 1063 do icell = 1,ncell 1064 ii = index_cells(4,icell); 1065 i1=index_cells(1,icell); i2=index_cells(2,icell); i3=index_cells(3,icell) 1066 ia1 = 0 ; ib1 = 0 1067 ! Loop over coefficients 1068 do icoeff=1,ncoeff 1069 ! Set the value of the coefficient 1070 coeff = coefficients(icoeff)%coefficient 1071 ! Loop over terms of this coefficient 1072 do iterm=1,coefficients(icoeff)%nterm 1073 ! Set the weight of this term 1074 weight =coefficients(icoeff)%terms(iterm)%weight 1075 tmp1 = one 1076 ndisp = coefficients(icoeff)%terms(iterm)%ndisp 1077 nstrain = coefficients(icoeff)%terms(iterm)%nstrain 1078 ndisp_tot = ndisp + nstrain 1079 ! Loop over displacement and strain 1080 do idisp1=1,ndisp_tot 1081 ! Set to one the acculation of forces and strain 1082 tmp2 = one 1083 tmp3 = zero 1084 1085 ! Strain case idisp > ndisp 1086 if (idisp1 > ndisp)then 1087 tmp3 = one 1088 ! Set the power_strain of the strain: 1089 idisp1_strain = idisp1 - ndisp 1090 power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp1_strain) 1091 ! Get the direction of the displacement or strain 1092 idir1 = coefficients(icoeff)%terms(iterm)%strain(idisp1_strain) 1093 if(abs(strain(idir1)) > tol10)then 1094 ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z)) 1095 tmp1 = tmp1 * (strain(idir1))**power_strain 1096 if(power_strain > 1) then 1097 ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...)) 1098 tmp3 = tmp3 * power_strain*(strain(idir1))**(power_strain-1) 1099 end if 1100 else 1101 tmp1 = zero 1102 if(power_strain > 1) then 1103 tmp3 = zero 1104 end if 1105 end if 1106 else 1107 ! Set the power_disp of the displacement: 1108 power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp1) 1109 ! Get the direction of the displacement or strain 1110 idir1 = coefficients(icoeff)%terms(iterm)%direction(idisp1) 1111 ! Displacement case idir = 1, 2 or 3 1112 ! indexes of the cell of the atom a 1113 cell_atoma1 = coefficients(icoeff)%terms(iterm)%cell(:,1,idisp1) 1114 if(cell_atoma1(1)/=0.or.cell_atoma1(2)/=0.or.cell_atoma1(3)/=0) then 1115 ! if the cell is not 0 0 0 we apply PBC: 1116 cell_atoma1(1) = i1 + cell_atoma1(1) 1117 cell_atoma1(2) = i2 + cell_atoma1(2) 1118 cell_atoma1(3) = i3 + cell_atoma1(3) 1119 call getPBCIndexes_supercell(cell_atoma1(1:3),sc_size(1:3)) 1120 ! index of the first atom (position in the supercell if the cell is not 0 0 0) 1121 ia1 = (cell_atoma1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 1122 & (cell_atoma1(2)-1)*sc_size(3)*natom_uc+& 1123 & (cell_atoma1(3)-1)*natom_uc+& 1124 & coefficients(icoeff)%terms(iterm)%atindx(1,idisp1) 1125 else 1126 ! index of the first atom (position in the supercell if the cell is 0 0 0) 1127 ia1 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp1) 1128 end if 1129 1130 ! indexes of the cell of the atom b (with PBC) same as ia1 1131 cell_atomb1 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp1) 1132 if(cell_atomb1(1)/=0.or.cell_atomb1(2)/=0.or.cell_atomb1(3)/=0) then 1133 cell_atomb1(1) = i1 + cell_atomb1(1) 1134 cell_atomb1(2) = i2 + cell_atomb1(2) 1135 cell_atomb1(3) = i3 + cell_atomb1(3) 1136 call getPBCIndexes_supercell(cell_atomb1(1:3),sc_size(1:3)) 1137 1138 ! index of the second atom in the (position in the supercell if the cell is not 0 0 0) 1139 ib1 = (cell_atomb1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 1140 & (cell_atomb1(2)-1)*sc_size(3)*natom_uc+& 1141 & (cell_atomb1(3)-1)*natom_uc+& 1142 & coefficients(icoeff)%terms(iterm)%atindx(2,idisp1) 1143 else 1144 ! index of the first atom (position in the supercell if the cell is 0 0 0) 1145 ib1 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp1) 1146 end if 1147 1148 ! Get the displacement for the both atoms 1149 disp1 = disp(idir1,ia1) 1150 disp2 = disp(idir1,ib1) 1151 1152 if(abs(disp1) > tol10 .or. abs(disp2)> tol10)then 1153 ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z)) 1154 tmp1 = tmp1 * (disp1-disp2)**power_disp 1155 if(power_disp > 1) then 1156 ! Accumulate forces for each displacement (\sum (Y(A_x-O_x)^Y-1(A_y-O_c)^Z+...)) 1157 tmp2 = tmp2 * power_disp*(disp1-disp2)**(power_disp-1) 1158 end if 1159 else 1160 tmp1 = zero 1161 if(power_disp > 1) then 1162 tmp2 = zero 1163 end if 1164 end if 1165 end if 1166 1167 do idisp2=1,ndisp_tot 1168 1169 if(idisp2 /= idisp1) then 1170 if (idisp2 > ndisp)then 1171 idisp2_strain = idisp2 - ndisp 1172 idir2 = coefficients(icoeff)%terms(iterm)%strain(idisp2_strain) 1173 ! Strain case 1174 ! Set the power_strain of the strain: 1175 power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp2_strain) 1176 ! Accumulate energy forces 1177 tmp2 = tmp2 * (strain(idir2))**power_strain 1178 ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...)) 1179 tmp3 = tmp3 * (strain(idir2))**power_strain 1180 else 1181 idir2 = coefficients(icoeff)%terms(iterm)%direction(idisp2) 1182 cell_atoma2=coefficients(icoeff)%terms(iterm)%cell(:,1,idisp2) 1183 if(cell_atoma2(1)/=0.or.cell_atoma2(2)/=0.or.cell_atoma2(3)/=0) then 1184 cell_atoma2(1) = i1 + cell_atoma2(1) 1185 cell_atoma2(2) = i2 + cell_atoma2(2) 1186 cell_atoma2(3) = i3 + cell_atoma2(3) 1187 call getPBCIndexes_supercell(cell_atoma2(1:3),sc_size(1:3)) 1188 ! index of the first atom (position in the supercell and direction) 1189 ! if the cell of the atom a is not 0 0 0 (may happen) 1190 ia2 = (cell_atoma2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 1191 & (cell_atoma2(2)-1)*sc_size(3)*natom_uc+& 1192 & (cell_atoma2(3)-1)*natom_uc+& 1193 & coefficients(icoeff)%terms(iterm)%atindx(1,idisp2) 1194 else 1195 ! index of the first atom (position in the supercell and direction) 1196 ia2 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp2) 1197 end if 1198 1199 cell_atomb2= coefficients(icoeff)%terms(iterm)%cell(:,2,idisp2) 1200 1201 if(cell_atomb2(1)/=0.or.cell_atomb2(2)/=0.or.cell_atomb2(3)/=0) then 1202 ! indexes of the cell2 (with PBC) 1203 cell_atomb2(1) = i1 + cell_atomb2(1) 1204 cell_atomb2(2) = i2 + cell_atomb2(2) 1205 cell_atomb2(3) = i3 + cell_atomb2(3) 1206 call getPBCIndexes_supercell(cell_atomb2(1:3),sc_size(1:3)) 1207 1208 ! index of the second atom in the (position in the supercell) 1209 ib2 = (cell_atomb2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 1210 & (cell_atomb2(2)-1)*sc_size(3)*natom_uc+& 1211 & (cell_atomb2(3)-1)*natom_uc+& 1212 & coefficients(icoeff)%terms(iterm)%atindx(2,idisp2) 1213 else 1214 ib2 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp2) 1215 end if 1216 1217 disp1 = disp(idir2,ia2) 1218 disp2 = disp(idir2,ib2) 1219 1220 ! Set the power_disp of the displacement: 1221 power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp2) 1222 tmp2 = tmp2 * (disp1-disp2)**power_disp 1223 tmp3 = tmp3 * (disp1-disp2)**power_disp 1224 1225 end if 1226 end if 1227 end do 1228 1229 if(idisp1 > ndisp)then 1230 ! Accumule stress tensor 1231 strten(idir1) = strten(idir1) + coeff * weight * tmp3 1232 else 1233 ! Accumule forces 1234 fcart(idir1,ia1) = fcart(idir1,ia1) + coeff * weight * tmp2 1235 fcart(idir1,ib1) = fcart(idir1,ib1) - coeff * weight * tmp2 1236 end if 1237 end do 1238 1239 energy_coeff(icoeff) = energy_coeff(icoeff) + coeff * weight * tmp1 1240 ! accumule energy 1241 energy = energy + coeff * weight * tmp1 1242 1243 end do 1244 end do 1245 end do 1246 1247 1248 ! MPI_SUM 1249 call xmpi_sum(energy, comm, ierr) 1250 call xmpi_sum(fcart , comm, ierr) 1251 call xmpi_sum(strten , comm, ierr) 1252 1253 !Write to anharmonic_energy_terms.out ORIGINAL 1254 INQUIRE(FILE=name_file,OPENED=file_opened,number=unit_out) 1255 if(file_opened .eqv. .TRUE.)then 1256 do icoeff=1,ncoeff 1257 call xmpi_sum(energy_coeff(icoeff), comm, ierr) 1258 ! Marcus write energy contributions of anharmonic terms to file 1259 if(icoeff <ncoeff)then 1260 write(unit_out,'(A,1ES24.16)',advance='no') ' ',energy_coeff(icoeff) 1261 else if(icoeff==ncoeff)then 1262 write(unit_out,'(A,1ES24.16)',advance='yes') ' ',energy_coeff(icoeff) 1263 end if 1264 enddo 1265 end if 1266 1267 1268 end subroutine polynomial_coeff_evaluate
m_polynomial_coeff/polynomial_coeff_free [ Functions ]
NAME
polynomial_coeff_free
FUNCTION
Free polynomial_coeff datatype
INPUTS
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
OUTPUT
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
SOURCE
308 subroutine polynomial_coeff_free(polynomial_coeff) 309 310 implicit none 311 312 !Arguments ------------------------------------ 313 !scalars 314 !arrays 315 type(polynomial_coeff_type), intent(inout) :: polynomial_coeff 316 !Local variables------------------------------- 317 !scalar 318 integer :: ii 319 !arrays 320 321 ! ************************************************************************* 322 323 if(allocated(polynomial_coeff%terms))then 324 do ii = 1,polynomial_coeff%nterm 325 call polynomial_term_free(polynomial_coeff%terms(ii)) 326 end do 327 ABI_FREE(polynomial_coeff%terms) 328 end if 329 polynomial_coeff%name = "" 330 polynomial_coeff%nterm = 0 331 polynomial_coeff%coefficient = zero 332 333 334 end subroutine polynomial_coeff_free
m_polynomial_coeff/polynomial_coeff_getEvenAnhaStrain [ Functions ]
NAME
polynomial_coeff_getEvenAnhaStrain
FUNCTION
Get even anharmonic strain terms in defined range of order
INPUTS
OUTPUT
polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with the polynomial_coeff ncoeff_out = number of coefficients
SOURCE
3934 subroutine polynomial_coeff_getEvenAnhaStrain(strain_terms,crystal,irred_ncoeff,power_strain,comm) 3935 3936 implicit none 3937 3938 !Arguments ------------------------------------ 3939 type(polynomial_coeff_type),allocatable,intent(inout) :: strain_terms(:) 3940 type(crystal_t), intent(inout) :: crystal 3941 integer,intent(out) :: irred_ncoeff 3942 integer,intent(in) :: power_strain(2) 3943 integer,intent(in) :: comm 3944 !scalars 3945 !arrays 3946 !Local variables------------------------------- 3947 real(dp) :: cutoff,coeff_ini 3948 integer :: ncoeff 3949 integer :: power_strph 3950 integer :: option 3951 integer :: icoeff1,icoeff2,start 3952 integer :: ncoeff_out 3953 logical,allocatable :: same(:) 3954 type(polynomial_coeff_type),allocatable :: strain_terms_tmp(:) 3955 !scalar 3956 !arrays 3957 integer :: sc_size(3) 3958 ! ************************************************************************* 3959 3960 !Initialize Variables for call to polynomial_coeff_getNorder 3961 cutoff = zero 3962 power_strph = zero 3963 option = 0 3964 sc_size = (/1,1,1/) 3965 coeff_ini = 1000000 3966 3967 3968 call polynomial_coeff_getNorder(strain_terms_tmp,crystal,cutoff,ncoeff,ncoeff_out,power_strain,& 3969 & power_strph,option,sc_size,comm,anharmstr=.true.,spcoupling=.false.,& 3970 & only_odd_power=.false.,only_even_power=.true.,compute_symmetric=.false.,& 3971 verbose=.false.) 3972 3973 3974 !TODO Probably put in one routine 3975 !Get irreducible strain 3976 ABI_MALLOC(same,(ncoeff_out)) 3977 same = .false. 3978 do icoeff1 =1,ncoeff_out 3979 start = icoeff1 + 1 3980 !write(*,*) "name coeff_",icoeff1,": ", strain_terms_tmp(icoeff1)%name 3981 do icoeff2=start,ncoeff_out 3982 if(.not.same(icoeff2)) same(icoeff2) = coeffs_compare(strain_terms_tmp(icoeff1),strain_terms_tmp(icoeff2)) 3983 enddo 3984 !write(*,*) "same(",icoeff1,"): ", same(icoeff1) 3985 enddo 3986 3987 irred_ncoeff = ncoeff_out - count(same) 3988 ABI_MALLOC(strain_terms,(irred_ncoeff)) 3989 3990 !Transfer irreducible strain to output array 3991 icoeff2=0 3992 icoeff1=0 3993 do icoeff1=1,ncoeff_out 3994 if(.not.same(icoeff1))then 3995 icoeff2=icoeff2 + 1 3996 call polynomial_coeff_init(coeff_ini,strain_terms_tmp(icoeff1)%nterm,strain_terms(icoeff2),& 3997 & strain_terms_tmp(icoeff1)%terms,strain_terms_tmp(icoeff1)%name,check=.TRUE.) 3998 endif 3999 enddo 4000 4001 !TEST MS write coefficients to xml to check 4002 !fname="EvenAnhaStrain.xml" 4003 !call polynomial_coeff_writeXML(strain_terms,irred_ncoeff,filename=fname,newfile=.true.) 4004 !write(*,*) "Strain terms to xml done" 4005 4006 4007 !Deallocateion 4008 call polynomial_coeff_list_free(strain_terms_tmp) 4009 ABI_FREE(same) 4010 4011 end subroutine polynomial_coeff_getEvenAnhaStrain
m_polynomial_coeff/polynomial_coeff_getList [ Functions ]
NAME
polynomial_coeff_getList
FUNCTION
Get the list of all the possible coefficients for the polynome
INPUTS
cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...) dist(3,natom,natom,nrpt) = distance between atoms atm1 is in the cell 0 0 0 atm2 is in the nrpt cell (see cell(3,nrpt)) for each component x,y and z crystal<type(crystal_t)> = datatype with all the information for the crystal natom = number of atoms in the unit cell nrpt = number of cell in the supercell
OUTPUT
list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients, for each coefficients (ncoeff_sym), we store the symmetrics(nsym) the 6th first dimensions are : 1 = direction of the IFC 2 = index of the atom number 1 (1=>natom) 3 = index of the atom number 2 (1=>natom) 4 = indexes of the cell of the second atom (the atom number 1 is always in the cell 0 0 0) 5 = weight of the term (-1 or 1) 6 = indexes of the symmetric list_symstr(nstr_sym,nsym) = array with the list of the strain and the symmetrics nstr_sym = number of coefficient for the strain ncoeff_sym = number of coefficient for the IFC range_ifc(3) = maximum cut-off for the inter atomic forces constants in each direction sc_size(3) = optional,size of the supercell used for the fit. For example if you want to fit 2x2x2 cell the interation Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process If check_pbc is true we remove these kind of terms
SOURCE
1423 subroutine polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,& 1424 & natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size,& 1425 & fit_iatom) 1426 1427 implicit none 1428 1429 !Arguments ------------------------------------ 1430 !scalars 1431 integer,intent(in) :: natom,nrpt 1432 real(dp), intent(in) :: cutoff 1433 integer,intent(out) :: ncoeff_sym,nstr_sym 1434 integer,optional,intent(in):: fit_iatom 1435 !arrays 1436 integer,intent(in) :: cell(3,nrpt) 1437 real(dp),intent(in):: dist(3,natom,natom,nrpt) 1438 type(crystal_t), intent(in) :: crystal 1439 integer,allocatable,intent(out) :: list_symcoeff(:,:,:),list_symstr(:,:,:) 1440 integer,optional,intent(in) :: sc_size(3) 1441 real(dp),intent(in):: range_ifc(3) 1442 !Local variables------------------------------- 1443 !scalar 1444 integer :: ia,ib,icoeff,icoeff2,icoeff_tot,icoeff_tmp,idisy1,idisy2,ii 1445 integer :: ipesy1,ipesy2,isym,irpt,irpt3,irpt_ref,irpt_sym 1446 integer :: jj,jsym,mu,fit_iatom_in 1447 integer :: ncoeff,ncoeff2,ncoeff3,ncoeff_max,nu 1448 integer :: nsym,shift_atm1(3) 1449 integer :: shift_atm2(3) 1450 real(dp):: dist_orig,dist_sym,tolsym8 1451 logical :: found,check_pbc,possible 1452 !arrays 1453 integer :: isym_rec(3,3),isym_rel(3,3),sc_size_in(3) 1454 integer :: transl(3),min_range(3),max_range(3) 1455 integer,allocatable :: blkval(:,:,:,:,:),list(:),list_symcoeff_tmp(:,:,:),list_symcoeff_tmp2(:,:,:) 1456 integer,allocatable :: list_symstr_tmp(:,:,:),indsym(:,:,:) ,symrec(:,:,:),symrel(:,:,:),list_symcoeff_tmp3(:,:,:) 1457 integer,allocatable :: index_irred(:) 1458 real(dp),allocatable :: tnons(:,:) 1459 real(dp),allocatable :: wkdist(:),distance(:,:,:) 1460 real(dp) :: difmin(3) 1461 real(dp) :: rprimd(3,3) 1462 real(dp) :: tratom(3) 1463 character(len=500) :: message 1464 1465 1466 1467 !Initialisation of variables 1468 irpt_sym = 0 1469 nsym = crystal%nsym 1470 rprimd = crystal%rprimd 1471 !ABI_MALLOC(xcart,(3,natom)) 1472 !ABI_MALLOC(xred,(3,natom)) 1473 !xcart(:,:) = crystal%xcart(:,:) 1474 !xred(:,:) = crystal%xred(:,:) 1475 ncoeff_max = nrpt*natom*natom*3*3 1476 1477 !Found the ref cell 1478 irpt_ref = 1 1479 do irpt=1,nrpt 1480 if(all(cell(:,irpt)==0))then 1481 irpt_ref = irpt 1482 ! exit 1483 end if 1484 end do 1485 1486 !Set the size of the interaction 1487 check_pbc = .FALSE. 1488 sc_size_in = 0 1489 min_range = 0; max_range = 0 1490 if(present(sc_size))then 1491 sc_size_in = sc_size 1492 do mu=1,3 1493 call findBound_supercell(min_range(mu),max_range(mu),sc_size_in(mu)) 1494 end do 1495 end if 1496 1497 !Check which atom to fit, if not present do all atoms 1498 if(present(fit_iatom))then 1499 fit_iatom_in = fit_iatom 1500 else 1501 fit_iatom_in = -1 1502 endif 1503 if(fit_iatom_in > natom)then 1504 write(message, '(3a)' )& 1505 & 'fit_iatom cannot be greater than the number of atoms on the reference unit cell',ch10,& 1506 & 'Action: Change input' 1507 ABI_ERROR(message) 1508 end if 1509 1510 !Obtain a list of rotated atom labels: 1511 ABI_MALLOC(indsym,(4,nsym,natom)) 1512 ABI_MALLOC(symrec,(3,3,nsym)) 1513 ABI_MALLOC(symrel,(3,3,nsym)) 1514 ABI_MALLOC(tnons,(3,nsym)) 1515 symrec = crystal%symrec 1516 symrel = crystal%symrel 1517 tnons = crystal%tnons 1518 1519 tolsym8=tol13 1520 call symatm(indsym,natom,nsym,symrec,tnons,& 1521 & tolsym8,crystal%typat,crystal%xred) 1522 ABI_MALLOC(blkval,(3,natom,3,natom,nrpt)) 1523 ABI_MALLOC(list,(natom*nrpt)) 1524 ABI_MALLOC(list_symcoeff_tmp,(5,ncoeff_max,nsym)) 1525 ABI_MALLOC(wkdist,(natom*nrpt)) 1526 1527 !1-Fill strain list 1528 ABI_MALLOC(list_symstr_tmp,(6,nsym,2)) 1529 list_symstr_tmp = 1 1530 do ia=1,6 1531 if(list_symstr_tmp(ia,1,1)==0)cycle 1532 ! Transform the voigt notation 1533 if(ia<=3)then 1534 mu=ia;nu=ia 1535 else 1536 select case(ia) 1537 case(4) 1538 mu=2;nu=3 1539 case(5) 1540 mu=1;nu=3 1541 case(6) 1542 mu=1;nu=2 1543 end select 1544 end if 1545 do isym=1,nsym 1546 ! Get the symmetry matrix 1547 isym_rel(:,:) = crystal%symrel(:,:,isym) 1548 do idisy1=1,3 1549 do idisy2=1,3 1550 if((isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)) then 1551 ! Transform to the voig notation 1552 if(idisy1==idisy2)then 1553 list_symstr_tmp(ia,isym,1) = idisy1 1554 list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1) 1555 else 1556 if(idisy1==1.or.idisy2==1)then 1557 if(idisy1==2.or.idisy2==2)then 1558 list_symstr_tmp(ia,isym,1) = 6 1559 end if 1560 if(idisy1==3.or.idisy2==3)then 1561 list_symstr_tmp(ia,isym,1) = 5 1562 end if 1563 else 1564 list_symstr_tmp(ia,isym,1) = 4 1565 end if 1566 end if 1567 list_symstr_tmp(ia,isym,2) = isym_rel(mu,idisy1) * isym_rel(nu,idisy2) 1568 end if 1569 end do 1570 end do 1571 ! Remove the symetric 1572 ! if(list_symstr_tmp(ia,isym,1) > ia) then 1573 ! list_symstr_tmp(list_symstr_tmp(ia,isym,1),:,1) = 0 1574 ! end if 1575 end do 1576 end do 1577 1578 !Count the number of strain and transfert into the final array 1579 nstr_sym = 0 1580 do ia=1,6 1581 if(list_symstr_tmp(ia,1,1)/=0) nstr_sym = nstr_sym + 1 1582 end do 1583 1584 if(allocated(list_symstr))then 1585 ABI_FREE(list_symstr) 1586 end if 1587 ABI_MALLOC(list_symstr,(nstr_sym,nsym,2)) 1588 1589 icoeff_tmp = 1 1590 do ia=1,6 1591 if(list_symstr_tmp(ia,1,1)/=0) then 1592 list_symstr(icoeff_tmp,:,:) = list_symstr_tmp(ia,:,:) 1593 icoeff_tmp = icoeff_tmp + 1 1594 end if 1595 end do 1596 !END STRAIN 1597 1598 !Compute the distance between each atoms. Indeed the dist array contains the difference of 1599 !cartesian coordinate for each direction 1600 ! dist: vector between atom a, (0,0,0) and atom b (rpt) 1601 1602 ABI_MALLOC(distance,(natom,natom,nrpt)) 1603 ! Fortran 2008: distance = norm2(dist, dim=1) 1604 do ia=1,natom 1605 do ib=1,natom 1606 do irpt=1,nrpt 1607 distance(ia,ib,irpt) = ((dist(1,ia,ib,irpt))**2+(dist(2,ia,ib,irpt))**2+& 1608 & (dist(3,ia,ib,irpt))**2)**0.5 1609 end do 1610 end do 1611 end do 1612 1613 1614 !Set to one blkval, all the coeff have to be compute 1615 blkval = 1 1616 icoeff = 1 1617 icoeff_tot = 1 1618 list_symcoeff_tmp = 0 1619 1620 !2-Fill atom list 1621 !Big loop over generic atom 1622 do ia=1,natom 1623 wkdist(:)=reshape(distance(ia,:,:),(/natom*nrpt/)) 1624 do ii=1,natom*nrpt 1625 list(ii)=ii 1626 end do 1627 ! FIXME: hexu: I think this should be improved. 1628 ! It seems that it depends on the specific order of the 1629 ! cell, and the implementation of sort algorithm. 1630 ! The sort should preserve the order if two distances are the same. 1631 ! In the future, the order of the cell might need to be unified 1632 ! Also it seems to depend on the cubic cell. 1633 call sort_dp(natom*nrpt,wkdist,list,tol8) 1634 do ii=1,natom*nrpt 1635 ! Get the irpt and ib 1636 irpt=(list(ii)-1)/natom+1 1637 ib=list(ii)-natom*(irpt-1) 1638 possible = .true. 1639 !Old way with the cut off 1640 if(cutoff < distance(ia,ib,irpt))then 1641 possible = .false. 1642 else 1643 ! in each direction jj, the d component d(jj) should be smaller than range_ifc(jj) 1644 do jj=1,3 1645 ! if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj) > tol10.or.& 1646 ! & abs(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj)) < tol10)then 1647 if(abs(dist(jj,ia,ib,irpt)) - range_ifc(jj) > tol10)then 1648 possible = .false. 1649 end if 1650 end do 1651 endif 1652 1653 ! If this distance is superior to the cutoff, we don't compute that term 1654 if(.not.possible)then 1655 blkval(:,ia,:,ib,irpt)= 0 1656 ! remove duplication: if in reference cell, only consider ab, not ba. 1657 if(irpt==irpt_ref)blkval(:,ib,:,ia,irpt)= 0 1658 ! Stop the loop 1659 cycle 1660 end if 1661 1662 ! If this coefficient is not possible, we cycle... 1663 if (all(blkval(:,ia,:,ib,irpt)==0)) cycle 1664 1665 ! Save the distance between the two atoms for futur checks 1666 dist_orig = (dist(1,ia,ib,irpt)**2+dist(2,ia,ib,irpt)**2+dist(3,ia,ib,irpt)**2)**0.5 1667 1668 do mu=1,3 1669 do nu=1,3 1670 ! Check if : - The coefficient is not yet compute 1671 ! - The directions are the same 1672 ! - The atoms are not equivalent 1673 if (mu/=nu) then 1674 blkval(mu,ia,nu,ib,irpt)=0 1675 blkval(nu,ia,mu,ib,irpt)=0 1676 cycle 1677 end if 1678 ! Pass if the atoms are identical and in the ref cell 1679 if(irpt==irpt_ref.and.ia==ib) then 1680 blkval(mu,ia,nu,ib,irpt)=0 1681 blkval(nu,ib,mu,ia,irpt)=0 1682 cycle 1683 end if 1684 if(blkval(mu,ia,nu,ib,irpt)==1)then 1685 ! Loop over symmetries 1686 do isym=1,nsym 1687 ! Get the symmetry matrix for this sym 1688 isym_rec(:,:) = crystal%symrec(:,:,isym) 1689 isym_rel(:,:) = crystal%symrel(:,:,isym) 1690 ! Get the corresponding atom and shift with the symetries 1691 ! For atom 1 1692 ipesy1 = indsym(4,isym,ia) 1693 shift_atm1 = indsym(1:3,isym,ia) 1694 ! And atom 2 1695 do jj=1,3 ! Apply transformation to original coordinates. 1696 tratom(jj) = dble(isym_rec(1,jj))*(crystal%xred(1,ib)+cell(1,irpt)-tnons(1,isym))& 1697 & +dble(isym_rec(2,jj))*(crystal%xred(2,ib)+cell(2,irpt)-tnons(2,isym))& 1698 & +dble(isym_rec(3,jj))*(crystal%xred(3,ib)+cell(3,irpt)-tnons(3,isym)) 1699 1700 end do 1701 1702 ! Find symmetrically equivalent atom 1703 call symchk(difmin,ipesy2,natom,tratom,transl,crystal%typat(ib),& 1704 & crystal%typat,crystal%xred(:,:)) 1705 1706 ! Put information into array indsym: translations and label 1707 shift_atm2(:)= transl(:) - shift_atm1(:) 1708 found = .false. 1709 do irpt3=1,nrpt 1710 if(cell(1,irpt3)==shift_atm2(1).and.& 1711 & cell(2,irpt3)==shift_atm2(2).and.& 1712 & cell(3,irpt3)==shift_atm2(3))then 1713 found = .true. 1714 irpt_sym = irpt3 1715 end if 1716 end do 1717 1718 ! Check the distance 1719 dist_sym = (dist(1,ipesy1,ipesy2,irpt_sym)**2+& 1720 & dist(2,ipesy1,ipesy2,irpt_sym)**2+& 1721 & dist(3,ipesy1,ipesy2,irpt_sym)**2)**0.5 1722 if(abs(dist_orig - dist_sym) > tol10)then 1723 write(message, '(a,i0,2a,I0,a,es15.8,2a,es15.8,2a)' )& 1724 & 'The distance between the atoms for the coefficient number ',icoeff,ch10,& 1725 & 'with the symmetry ',isym,' is ',dist_sym,ch10,'but the original distance is',& 1726 & dist_orig,ch10,& 1727 & 'Action: Contact abinit group' 1728 ABI_BUG(message) 1729 end if 1730 ! Now that a symmetric perturbation has been obtained, 1731 ! including the expression of the symmetry matrix, see 1732 ! if the symmetric perturbations are available 1733 do idisy1=1,3 1734 do idisy2=1,3 1735 if (idisy1/=idisy2) then 1736 ! Remove this term (is not computed) 1737 ! Also remove opposite term... (Srx-Tix) = (Tix-Srx) 1738 blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0 1739 blkval(idisy2,ipesy1,idisy1,ipesy2,irpt_sym) = 0 1740 cycle 1741 else 1742 if(isym_rel(mu,idisy1)/=0.and.isym_rel(nu,idisy2)/=0)then 1743 if(.not.found.or.(irpt_sym==irpt_ref.and.ipesy1==ipesy2)) then 1744 ! Remove this term (is not computed) Sr-Sr or not include in the cell 1745 ! Also remove oposite term... (Srx-Tix) = (Ti-Srx) 1746 blkval(idisy1,ipesy1,idisy2,ipesy2,irpt_sym) = 0 1747 blkval(idisy2,ipesy2,idisy1,ipesy1,irpt_sym) = 0 1748 cycle 1749 else 1750 ! Fill the list with the coeff and symmetric (need all symmetrics) 1751 list_symcoeff_tmp(1:4,icoeff,isym)=(/idisy1,ipesy1,ipesy2,irpt_sym/) 1752 ! Check the sign 1753 if(isym_rel(mu,idisy1)/=isym_rel(nu,idisy2))then 1754 write(message, '(a,i0,a,I0,4a)' )& 1755 & 'The sign of coefficient number ',icoeff,' with the symmetry ',isym,ch10,& 1756 & 'can not be found... Something is going wrong',ch10,& 1757 & 'Action: Contact abinit group' 1758 ABI_BUG(message) 1759 end if 1760 list_symcoeff_tmp(5,icoeff,isym)= isym_rel(nu,idisy2) 1761 end if 1762 end if 1763 end if 1764 end do 1765 end do 1766 end do ! end loop sym 1767 icoeff = icoeff + 1 1768 end if 1769 ! This coeff is now computed 1770 blkval(mu,ia,nu,ib,irpt)= 0 1771 end do ! end loop nu 1772 end do ! end loop mu 1773 end do ! end loop ii: iatom and irpt. 1774 end do ! end loop ia 1775 1776 !Reset the output 1777 if(allocated(list_symcoeff))then 1778 ABI_FREE(list_symcoeff) 1779 end if 1780 1781 ABI_FREE(distance) 1782 1783 !Transfert the final array with all the coefficients 1784 !With this array, we can access to all the terms presents 1785 !ncoeff1 + symetrics 1786 !first dimension is 4 (mu,ia,ib,irpt) 1787 !irpt is the index of the cell of the atom ib in the cell array 1788 !example cell(:,irpt=12) can be (-1 0 -2). The cell of ia is 1789 !always 0 0 0 1790 !Transfert the final array for the list of irreductible coeff and symetries 1791 !With this array, we can access to the irretuctible coefficients (ncoeff1) and 1792 !all the symetrics of these coefficients (nsym) 1793 !first dimension is 5 (mu,ia,ib,irpt,icoeff) 1794 !icoeff is the position of this coefficients in the list_fullcoeff array 1795 1796 !1/ step remove the zero coeff in this array 1797 ncoeff = 0 1798 do icoeff = 1,ncoeff_max 1799 if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then 1800 ncoeff = ncoeff + 1 1801 end if 1802 end do 1803 1804 ABI_MALLOC(list_symcoeff_tmp2,(6,ncoeff,nsym)) 1805 list_symcoeff_tmp2 = 0 1806 icoeff = 0 1807 do icoeff_tmp = 1,ncoeff_max 1808 if(.not.(all(list_symcoeff_tmp(:,icoeff_tmp,1)==0)))then 1809 icoeff = icoeff + 1 1810 list_symcoeff_tmp2(1:5,icoeff,:) = list_symcoeff_tmp(1:5,icoeff_tmp,:) 1811 end if 1812 end do 1813 1814 1815 !2/ set the dimension six of list_symcoeff_tmp2(6,icoeffs,1) 1816 ! and check is a symetric coeff is not coresspondig to an other 1817 ! one, in this case we set this coeff to 0 1818 ! ncoeff2 = zero 1819 do icoeff = 1,ncoeff 1820 ! found the index of each coeff in list_fullcoeff 1821 do isym = 1,nsym 1822 icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,1),& 1823 & list_symcoeff_tmp2(2,icoeff,isym),& 1824 & list_symcoeff_tmp2(3,icoeff,isym),& 1825 & list_symcoeff_tmp2(4,icoeff,isym),& 1826 & list_symcoeff_tmp2(1,icoeff,isym),& 1827 & ncoeff) 1828 list_symcoeff_tmp2(6,icoeff,isym) = icoeff2 1829 end do 1830 end do 1831 1832 !2.5/do checks 1833 do icoeff = 1,ncoeff 1834 do isym = 1,nsym 1835 if(list_symcoeff_tmp2(6,icoeff,isym)==0)then 1836 write(message, '(a,i0,a,I0,4a)' )& 1837 & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,& 1838 & 'has no equivalent',ch10,& 1839 & 'Action: Contact abinit group' 1840 ABI_BUG(message) 1841 else 1842 if(icoeff /= list_symcoeff_tmp2(6,icoeff,isym))then 1843 if(list_symcoeff_tmp2(1,icoeff,isym)/=& 1844 & list_symcoeff_tmp2(1,list_symcoeff_tmp2(6,icoeff,isym),1))then 1845 write(message, '(a,i0,a,I0,2a,I0,4a)' )& 1846 & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,& 1847 & 'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,& 1848 & 'because the direction is different:',ch10,& 1849 & 'Action: Contact abinit group' 1850 ABI_BUG(message) 1851 end if 1852 if(list_symcoeff_tmp2(4,icoeff,isym)/=& 1853 & list_symcoeff_tmp2(4,list_symcoeff_tmp2(6,icoeff,isym),1))then 1854 write(message, '(a,i0,a,I0,2a,I0,4a)' )& 1855 & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,& 1856 & 'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,& 1857 & 'because the cell is different',ch10,& 1858 & 'Action: Contact abinit group' 1859 ABI_BUG(message) 1860 end if 1861 if((list_symcoeff_tmp2(2,icoeff,isym)/=& 1862 & list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1).and.& 1863 & list_symcoeff_tmp2(3,icoeff,isym)/=& 1864 & list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1)).and.& 1865 & (list_symcoeff_tmp2(2,icoeff,isym)/=& 1866 & list_symcoeff_tmp2(3,list_symcoeff_tmp2(6,icoeff,isym),1).and.& 1867 & list_symcoeff_tmp2(3,icoeff,isym)/=& 1868 & list_symcoeff_tmp2(2,list_symcoeff_tmp2(6,icoeff,isym),1)))then 1869 write(message, '(a,i0,a,I0,2a,I0,4a)' )& 1870 & 'The coefficient number ',icoeff,' with the symetrie ',isym,ch10,& 1871 & 'does not refer to the same coefficient ',list_symcoeff_tmp2(6,icoeff,1),ch10,& 1872 & 'because the atoms different',ch10,& 1873 & 'Action: Contact abinit group' 1874 ABI_BUG(message) 1875 end if 1876 end if 1877 end if 1878 end do 1879 end do 1880 1881 1882 !Check if the atom 2 is not (with the PBC) is in the same cell than 1883 !the atom 1. For example if you want to fit 2x2x2 cell the interation 1884 !Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process... 1885 do icoeff = 1,ncoeff 1886 if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then 1887 do mu=1,3 1888 ! out of max_range(mu) 1889 if( -max_range(mu) > cell(mu,list_symcoeff_tmp2(4,icoeff,1)) .or. & 1890 & cell(mu,list_symcoeff_tmp2(4,icoeff,1)) > max_range(mu))then 1891 1892 ! FIXME: hexu: cell(mu, irpt). But why -1? 1893 ! The if above looks sufficient. 1894 !if(cell(mu,list_symcoeff_tmp2(4,icoeff,1)) < -1)then 1895 list_symcoeff_tmp2(:,icoeff,:)=0 1896 exit 1897 !end if 1898 end if 1899 end do 1900 end if 1901 !MS only keep terms with ia == fit_iatom or ib == fit_iatom specified in input 1902 if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)) .and. fit_iatom_in > 0)then 1903 if (list_symcoeff_tmp2(2,icoeff,1) /= fit_iatom_in) then !& LB 1904 !& list_symcoeff_tmp2(3,icoeff,1) /= fit_iatom_in)then LB 1905 list_symcoeff_tmp2(:,icoeff,1) = 0 1906 ! else if(list_symcoeff_tmp2(2,icoeff,1) /= fit_iatom_in) then ! & 1907 !& .and. list_symcoeff_tmp2(3,icoeff,1) == fit_iatom_in & !LB 1908 !& .and. any(cell(:,list_symcoeff_tmp2(4,icoeff,1)) /= 0))then !LB 1909 ! list_symcoeff_tmp2(:,icoeff,1) = 0 1910 endif 1911 endif 1912 !For Debugging keep only terms A_x-A_x[100] etc. comment if above 1913 ! if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then 1914 ! if (list_symcoeff_tmp2(2,icoeff,1) /= 1 .or. list_symcoeff_tmp2(3,icoeff,1) /= 1)then 1915 ! list_symcoeff_tmp2(:,icoeff,1) = 0 1916 ! endif !MS only keep terms with ia == 1 1917 ! endif 1918 end do 1919 1920 !write(std_out,*) "DEBUG: min_range(:): ", min_range(:) 1921 !write(std_out,*) "DEBUG: max_range(:): ", max_range(:) 1922 1923 !3/ Remove useless terms like opposites 1924 do icoeff = 1,ncoeff 1925 if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)) .and. & ! valid term. 1926 & list_symcoeff_tmp2(2,icoeff,1) /= list_symcoeff_tmp2(3,icoeff,1))then ! on-site term. ! FIXME: hexu:I don't understand if iatom=jatom but rpt/=0, it is not onsite. 1927 do isym = 1,nsym 1928 !icoeff2 = list_symcoeff_tmp2(6,icoeff,isym) 1929 !if (icoeff2> icoeff)then 1930 ! list_symcoeff_tmp2(:,icoeff2,1) = 0 1931 !end if 1932 do jsym=1,nsym 1933 ! FIXME: Should the rpt be -rpt?? 1934 ! icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,jsym),& 1935 !& list_symcoeff_tmp2(3,icoeff,isym),& 1936 !& list_symcoeff_tmp2(2,icoeff,isym),& 1937 !& list_symcoeff_tmp2(4,icoeff,isym),& !rpt 1938 !& list_symcoeff_tmp2(1,icoeff,isym),& 1939 !& ncoeff) 1940 icoeff2 = getCoeffFromList(list_symcoeff_tmp2(:,:,jsym),& 1941 & list_symcoeff_tmp2(3,icoeff,isym),& 1942 & list_symcoeff_tmp2(2,icoeff,isym),& 1943 & list_symcoeff_tmp2(4,icoeff,isym),& !rpt 1944 & find_opposite_irpt(cell, list_symcoeff_tmp2(1,icoeff,isym)),& 1945 & ncoeff) 1946 1947 if (icoeff2> icoeff)then 1948 list_symcoeff_tmp2(:,icoeff2,1) = 0 1949 end if 1950 end do 1951 end do 1952 end if 1953 end do 1954 1955 !4/ Recount the number of coeff after step 3 1956 ncoeff2 = 0 1957 isym = 0 1958 do icoeff = 1,ncoeff 1959 if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0)))then 1960 ncoeff2 = ncoeff2 + 1 1961 end if 1962 end do 1963 1964 !write(std_out,*) "DEBUG ncoeff after 2.5.3: ", ncoeff 1965 !write(std_out,*) "DEBUG ncoeff2 after 2.5.3: ", ncoeff2 1966 1967 ABI_FREE(list_symcoeff_tmp) 1968 ABI_MALLOC(list_symcoeff_tmp,(6,ncoeff,nsym)) 1969 list_symcoeff_tmp = list_symcoeff_tmp2 1970 1971 !4.1 Count irreducible terms 1972 do icoeff = 1,ncoeff 1973 do isym = 1,nsym 1974 icoeff2 = list_symcoeff_tmp2(6,icoeff,isym) 1975 if (icoeff2> icoeff)then 1976 list_symcoeff_tmp(:,icoeff2,1) = 0 1977 end if 1978 end do 1979 end do 1980 1981 ncoeff3 = 0 1982 do icoeff = 1,ncoeff 1983 if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then 1984 ncoeff3 = ncoeff3 + 1 1985 end if 1986 end do 1987 !write(std_out,*) "DEBUG ncoeff3 after 2.5.3: ", ncoeff3 1988 1989 !4.2 Put irreducible terms in front of list_symcoeff_tmp3 1990 ! Store index of irreducible terms 1991 ABI_MALLOC(list_symcoeff_tmp3,(6,ncoeff2,nsym)) 1992 ABI_MALLOC(index_irred,(ncoeff3)) 1993 icoeff_tmp = 0 1994 do icoeff = 1,ncoeff 1995 if(.not.(all(list_symcoeff_tmp(:,icoeff,1)==0)))then 1996 icoeff_tmp = icoeff_tmp+1 1997 index_irred(icoeff_tmp) = icoeff 1998 list_symcoeff_tmp3(:,icoeff_tmp,:) = list_symcoeff_tmp(:,icoeff,:) 1999 endif 2000 enddo 2001 2002 !4.3 Put symmetric equivalents behind in list_symcoeff_tmp3 2003 ! Attention icoeff_tmps keeps it's value of loop before 2004 ! TODO for check should be equal to ncoeff2 2005 do icoeff = 1,ncoeff 2006 if(.not.(all(list_symcoeff_tmp2(:,icoeff,1)==0))& 2007 & .and..not. any(index_irred == icoeff))then 2008 icoeff_tmp = icoeff_tmp+1 2009 list_symcoeff_tmp3(:,icoeff_tmp,:) = list_symcoeff_tmp2(:,icoeff,:) 2010 endif 2011 enddo 2012 2013 !4.4 A little copy round 2014 ABI_FREE(list_symcoeff_tmp) 2015 ABI_MALLOC(list_symcoeff_tmp,(6,ncoeff2,nsym)) 2016 list_symcoeff_tmp = list_symcoeff_tmp3 2017 2018 !5/ Final transfert 2019 ABI_MALLOC(list_symcoeff,(6,ncoeff2,nsym)) 2020 list_symcoeff = 0 2021 icoeff = 0 2022 do icoeff = 1,ncoeff2 2023 list_symcoeff(1:6,icoeff,:) = list_symcoeff_tmp(1:6,icoeff,:) 2024 do isym=1,nsym 2025 end do 2026 end do 2027 2028 !6/ reset the dimension six of list_symcoeff_tmp2(6,icoeffs,1) 2029 ! and check is a symetric coeff is not coresspondig to an other 2030 ! one, in this case we set this coeff to 0 2031 do icoeff = 1,ncoeff2 2032 ! found the index of each coeff in list_fullcoeff 2033 do isym = 1,nsym 2034 icoeff2 = getCoeffFromList(list_symcoeff(:,:,1),& 2035 & list_symcoeff(2,icoeff,isym),& 2036 & list_symcoeff(3,icoeff,isym),& 2037 & list_symcoeff(4,icoeff,isym),& 2038 & list_symcoeff(1,icoeff,isym),& 2039 & ncoeff2) 2040 list_symcoeff(6,icoeff,isym) = icoeff2 2041 ! list_symcoeff(6,icoeff2,isym) = icoeff 2042 end do 2043 end do 2044 2045 !Set the max number of coeff inside list_symcoeff 2046 ncoeff_sym = ncoeff3 2047 2048 2049 !Deallocation 2050 ABI_FREE(blkval) 2051 ABI_FREE(list) 2052 ABI_FREE(list_symcoeff_tmp) 2053 ABI_FREE(list_symcoeff_tmp2) 2054 ABI_FREE(list_symcoeff_tmp3) 2055 ABI_FREE(list_symstr_tmp) 2056 ABI_FREE(index_irred) 2057 ABI_FREE(indsym) 2058 ABI_FREE(symrec) 2059 ABI_FREE(symrel) 2060 ABI_FREE(tnons) 2061 !ABI_FREE(xcart) 2062 !ABI_FREE(xred ) 2063 ABI_FREE(wkdist) 2064 end subroutine polynomial_coeff_getList
m_polynomial_coeff/polynomial_coeff_getName [ Functions ]
NAME
polynomial_coeff_getName
FUNCTION
get the name of a polynomial coefficient
INPUTS
natom = number of atoms polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype symbols(natom) = array with the atomic symbol:["Sr","Ru","O1","O2","O3"] recompute = (optional) flag to set if the name has to be recomputed iterm = (optional) number of the term used for the name
OUTPUT
name = name xof the coefficients
SOURCE
473 subroutine polynomial_coeff_getName(name,polynomial_coeff,symbols,recompute,iterm) 474 475 implicit none 476 477 !Arguments ------------------------------------ 478 !scalars 479 integer,optional,intent(in) :: iterm 480 !arrays 481 character(len=5),intent(in) :: symbols(:) 482 character(len=200),intent(out):: name 483 type(polynomial_coeff_type),optional, intent(in) :: polynomial_coeff 484 logical,optional,intent(in) :: recompute 485 !Local variables------------------------------- 486 !scalar 487 integer :: ii,idisp,iterm_in 488 logical :: need_recompute 489 !arrays 490 integer :: cell_atm1(3),cell_atm2(3) 491 character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/) 492 character(len=1) :: dir 493 character(len=2) :: power_disp,power_dispchar 494 character(len=20) :: atm1,atm2 495 character(len=100):: atm1_tmp,atm2_tmp 496 character(len=200):: text 497 character(len=500):: message 498 ! ************************************************************************* 499 500 !Reset output 501 name="" 502 iterm_in = 1 503 504 !Set the optional arguments 505 need_recompute = .FALSE. 506 if(present(recompute)) need_recompute = recompute 507 if(present(iterm)) then 508 iterm_in = iterm 509 else 510 if(need_recompute)then 511 iterm_in = -1 512 do ii=1,polynomial_coeff%nterm 513 ! Find the index of the ref 514 if(iterm_in==-1) then !Need to find the reference term 515 do idisp=1,polynomial_coeff%terms(ii)%ndisp 516 if(polynomial_coeff%terms(ii)%direction(idisp) > 0) then 517 iterm_in = ii 518 if(any(polynomial_coeff%terms(ii)%cell(:,1,idisp) /= 0).or.& 519 & any(polynomial_coeff%terms(ii)%cell(:,2,idisp) /= 0)) then 520 iterm_in = -1 521 exit 522 end if 523 end if 524 end do!end do disp 525 else 526 exit 527 end if 528 end do!end do term 529 ! If not find, we set to the first element 530 if(iterm_in==-1) iterm_in = 1 531 else 532 iterm_in = 1 533 end if 534 end if 535 !Do check 536 if(iterm_in > polynomial_coeff%nterm.or.iterm_in < 0) then 537 write(message, '(5a)')& 538 & ' The number of the requested term for the generation of',ch10,& 539 & 'the name of the coefficient is not possible.',ch10,& 540 & 'Action: Contact Abinit group.' 541 ABI_BUG(message) 542 end if 543 544 if(polynomial_coeff%name /= "".and..not.need_recompute)then 545 name = polynomial_coeff%name 546 else 547 ! Nedd to recompute 548 do idisp=1,polynomial_coeff%terms(iterm_in)%ndisp 549 text = "" 550 !Fill variables for this displacement 551 write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_disp(idisp) 552 power_disp=trim(power_dispchar) 553 554 atm1=symbols(polynomial_coeff%terms(iterm_in)%atindx(1,idisp)) 555 atm2=symbols(polynomial_coeff%terms(iterm_in)%atindx(2,idisp)) 556 dir=mutodir(polynomial_coeff%terms(iterm_in)%direction(idisp)) 557 cell_atm1=polynomial_coeff%terms(iterm_in)%cell(:,1,idisp) 558 cell_atm2=polynomial_coeff%terms(iterm_in)%cell(:,2,idisp) 559 ! Construct ATM1 560 if (any(cell_atm1(:) /= 0) )then 561 write(atm1_tmp,'(4a,I0,a,I0,a,I0,a)') trim(atm1),"_",dir,"[",cell_atm1(1)," ",& 562 & cell_atm1(2)," ",cell_atm1(3),"]" 563 else 564 atm1_tmp = trim(atm1)//"_"//dir 565 end if 566 ! Construct ATM2 567 if(any(cell_atm2(:) /= 0))then 568 write(atm2_tmp,'(4a,I0,a,I0,a,I0,a)') trim(atm2),"_",dir,"[",cell_atm2(1)," ",& 569 & cell_atm2(2)," ",cell_atm2(3),"]" 570 else 571 atm2_tmp = trim(atm2)//"_"//dir 572 end if 573 text="("//trim(atm1_tmp)//"-"//trim(atm2_tmp)//")^"//power_disp 574 name = trim(name)//trim(text) 575 end do 576 !Strain case 577 do idisp=1,polynomial_coeff%terms(iterm_in)%nstrain 578 write(power_dispchar,'(I0)') polynomial_coeff%terms(iterm_in)%power_strain(idisp) 579 power_disp=trim(power_dispchar) 580 dir=mutodir(3+polynomial_coeff%terms(iterm_in)%strain(idisp)) 581 text="("//"eta_"//trim(dir)//")^"//power_disp 582 name = trim(name)//trim(text) 583 end do 584 end if 585 586 end subroutine polynomial_coeff_getName
m_polynomial_coeff/polynomial_coeff_getNorder [ Functions ]
NAME
polynomial_coeff_getNorder
FUNCTION
Compute and store into the datatype coefficients, all the possible coefficients for given orders
INPUTS
cutoff = cut-off for the inter atomic forces constants crystal<type(crystal_t)> = datatype with all the information for the crystal power_disps(2) = array with the minimal and maximal power_disp to be computed max_power_strain = maximum order of the strain of the strain phonon coupling option = 0 compute all terms 1 still in development sc_size(3) = size of the supercell used for the fit. For example if you want to fit 2x2x2 cell the interation Sr-Ti and Sr-Ti[2 0 0] will be identical for the fit process If check_pbc is true we remove these kind of terms comm = MPI communicator anharmstr = logical, optional : TRUE, the anharmonic strain is computed (\eta)^power_disp ... FALSE, (default) the anharmonic strain are not computed spcoupling= logical, optional : TRUE(default) the anharmonic strain-phonon coupling is computed only_odd_power = logical, optional : if TRUE return only odd power only_even_power= logical, optional : if TRUe return only even power distributed = logical, optional : True, the coefficients will be distributed on the CPU verbose = optional, flag for the verbose mode
OUTPUT
polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff) = array of datatype with the polynomial_coeff ncoeff = number of coefficients for this CPU if distributed == true, all otherwise ncoeff_tot = total number of coefficient over the CPU
SOURCE
2104 subroutine polynomial_coeff_getNorder(coefficients,crystal,cutoff,ncoeff,ncoeff_tot,power_disps,& 2105 & max_power_strain,option,sc_size,comm,anharmstr,spcoupling,& 2106 & distributed,only_odd_power,only_even_power,fit_iatom,& 2107 & compute_symmetric,dispterms,verbose) 2108 2109 implicit none 2110 2111 !Arguments ------------------------------------ 2112 !scalars 2113 integer,intent(in) :: max_power_strain,option,comm 2114 integer,intent(out):: ncoeff,ncoeff_tot 2115 integer,optional,intent(in) :: fit_iatom 2116 real(dp),intent(in):: cutoff 2117 logical,optional,intent(in) :: anharmstr,spcoupling,distributed,verbose,dispterms 2118 logical,optional,intent(in) :: only_odd_power,only_even_power,compute_symmetric 2119 !arrays 2120 integer,intent(in) :: power_disps(2),sc_size(3) 2121 type(crystal_t), intent(inout) :: crystal 2122 type(polynomial_coeff_type),allocatable,intent(inout) :: coefficients(:) 2123 !Local variables------------------------------- 2124 !scalar 2125 integer :: icoeff,icoeff2,icoeff3,ierr,ii,iterm 2126 integer :: i 2127 integer :: master,my_rank,my_ncoeff,my_newncoeff,natom,ncombination,ncoeff_max,ncoeff_sym 2128 integer :: ncoeff_symsym,nirred_comb,iirred_comb,ndisp,nstrain,fit_iatom_in 2129 integer :: ncoeff_alone,ndisp_max,nproc,nrpt,nsym,nterm,nstr_sym,my_size 2130 integer :: my_icoeff,rank_to_send,rank_to_receive,rank_to_send_save 2131 integer :: ncombi_alone,my_ncombi_simple,my_ncombi_start,my_ncombi_end,my_ncombi,my_nirred 2132 logical :: iam_master,need_anharmstr,need_spcoupling,need_distributed,need_verbose 2133 logical :: need_only_odd_power,need_only_even_power,compute_sym,irreducible,need_compute_symmetric 2134 logical :: need_dispterms 2135 !arrays 2136 integer :: shape_listsymcoeff(3),shape_listsymstr(3) 2137 integer,allocatable :: buffsize(:),buffdispl(:) !,dummylist(:) ,index_irred(:) 2138 integer,allocatable :: offsets(:) 2139 integer,allocatable :: cell(:,:),compatibleCoeffs(:,:) 2140 integer,allocatable :: list_symcoeff(:,:,:),list_symstr(:,:,:),list_coeff(:),list_combination(:,:) 2141 integer,allocatable :: list_combination_tmp(:,:) 2142 integer,allocatable :: irank_ncombi(:),my_index_irredcomb(:) 2143 integer,allocatable :: my_coefflist(:),my_coeffindexes(:),my_newcoeffindexes(:),my_list_combination(:,:) 2144 type(int2d_array_type) :: my_array_combination 2145 integer,allocatable :: my_list_combination_tmp(:,:) 2146 real(dp) :: rprimd(3,3),range_ifc(3) 2147 real(dp),allocatable :: dist(:,:,:,:) 2148 character(len=5),allocatable :: symbols(:) 2149 character(len=200):: name 2150 character(len=500) :: message 2151 type(polynomial_coeff_type),dimension(:),allocatable :: coeffs_tmp 2152 type(polynomial_term_type),dimension(:),allocatable :: terms 2153 character(len=fnlen) :: filename 2154 ! ************************************************************************* 2155 !Hide filename for debugging 2156 ABI_UNUSED(filename) 2157 2158 !MPI variables 2159 master = 0 2160 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2161 iam_master = (my_rank == master) 2162 2163 !Free the output 2164 if(allocated(coefficients))then 2165 do ii =1,size(coefficients) 2166 call polynomial_coeff_free(coefficients(ii)) 2167 end do 2168 ABI_FREE(coefficients) 2169 end if 2170 2171 !Check 2172 if(option > power_disps(2))then 2173 write(message, '(3a)' )& 2174 & 'Option can not be superior to the maximum order ',ch10,& 2175 & 'Action: contact abinit group' 2176 ABI_ERROR(message) 2177 end if 2178 2179 !Initialisation of variables 2180 need_anharmstr = .TRUE. 2181 if(present(anharmstr)) need_anharmstr = anharmstr 2182 need_spcoupling = .TRUE. 2183 if(present(spcoupling)) need_spcoupling = spcoupling 2184 need_distributed = .FALSE. 2185 if(present(distributed)) need_distributed = distributed 2186 need_verbose = .TRUE. 2187 if(present(verbose)) need_verbose = verbose 2188 need_only_odd_power = .FALSE. 2189 if(present(only_odd_power)) need_only_odd_power = only_odd_power 2190 need_only_even_power = .FALSE. 2191 if(present(only_even_power)) need_only_even_power = only_even_power 2192 need_compute_symmetric = .TRUE. 2193 if(present(compute_symmetric)) need_compute_symmetric = compute_symmetric 2194 need_dispterms = .TRUE. 2195 if(present(dispterms)) need_dispterms = dispterms 2196 2197 if(need_only_odd_power.and.need_only_even_power)then 2198 write(message, '(3a)' )& 2199 & 'need_only_odd_power and need_only_even_power are both true',ch10,& 2200 & 'Action: contact abinit group' 2201 ABI_ERROR(message) 2202 end if 2203 2204 !Check which atom to fit, if not present do all atoms 2205 if(present(fit_iatom))then 2206 ! FIXME: the fit_iatom as the input is the index of the atom 2207 ! but in the fitting subroutine, it is the irreducible atom. 2208 fit_iatom_in = fit_iatom 2209 else 2210 fit_iatom_in = -1 2211 endif 2212 2213 rprimd = crystal%rprimd 2214 2215 != 2216 call prepare_for_getList(crystal,sc_size, dist, cell, natom, nsym, nrpt, range_ifc , symbols) 2217 2218 if(iam_master)then 2219 if(need_verbose)then 2220 write(message,'(1a)')' Generation of the list of all the possible pairs of atoms within cutoff' 2221 call wrtout(std_out,message,'COLL') 2222 end if 2223 call polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,& 2224 & natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size=sc_size,& 2225 & fit_iatom=fit_iatom_in) 2226 2227 2228 shape_listsymcoeff = shape(list_symcoeff) 2229 shape_listsymstr = shape(list_symstr) 2230 endif!if iam master 2231 2232 !Broadcast Results of getList 2233 call xmpi_bcast(shape_listsymcoeff, master, comm, ierr) 2234 call xmpi_bcast(shape_listsymstr, master, comm, ierr) 2235 call xmpi_bcast(nstr_sym, master, comm, ierr) 2236 call xmpi_bcast(ncoeff_sym, master, comm, ierr) 2237 if(.not. iam_master )then 2238 ABI_MALLOC(list_symcoeff,(shape_listsymcoeff(1),shape_listsymcoeff(2),shape_listsymcoeff(3))) 2239 ABI_MALLOC(list_symstr,(shape_listsymstr(1),shape_listsymstr(2),shape_listsymstr(3))) 2240 endif 2241 call xmpi_bcast(list_symcoeff, master, comm, ierr) 2242 call xmpi_bcast(list_symstr, master, comm, ierr) 2243 call xmpi_barrier(comm) 2244 2245 ncoeff_symsym = size(list_symcoeff(1,:,1)) 2246 2247 !Compute the total number of coefficient 2248 ncoeff_tot = ncoeff_sym+nstr_sym 2249 2250 !if(iam_master)then 2251 !Check the distanceance bewteen coefficients and store integer: 2252 ! 0: the mix between these coefficient is not possible 2253 ! 1: the mix between these coefficient is possible 2254 ABI_MALLOC(compatibleCoeffs,(ncoeff_symsym+nstr_sym,ncoeff_symsym+nstr_sym)) 2255 compatibleCoeffs(:,:) = 1 2256 2257 if(need_verbose)then 2258 write(message,'(1a)')' Check the compatible coefficients with respect to the cutoff' 2259 call wrtout(std_out,message,'COLL') 2260 end if 2261 2262 do icoeff=1,ncoeff_symsym+nstr_sym 2263 do icoeff2=1,ncoeff_symsym+nstr_sym 2264 ! Select case: 2265 ! if both icoeff are displacement => check the distance 2266 ! if both icoeff are strain => check the flag 2267 ! Otherwise cycle (we keep the term) 2268 ! if(icoeff>ncoeff_sym.and.icoeff2<=ncoeff_sym)cycle 2269 ! if(icoeff<=ncoeff_sym.and.icoeff2>ncoeff_sym)cycle 2270 !Forbid anharmonic strain terms if not wanted 2271 if((icoeff>ncoeff_symsym .and. icoeff2>ncoeff_symsym).and.& 2272 & .not.need_anharmstr) then 2273 compatibleCoeffs(icoeff,icoeff2) = 0 2274 compatibleCoeffs(icoeff2,icoeff) = 0 2275 end if 2276 !Forbid strain phonon terms if not wanted 2277 if((icoeff>ncoeff_symsym.or.icoeff2>ncoeff_symsym).and.& 2278 & .not.need_spcoupling) then 2279 compatibleCoeffs(icoeff,icoeff2) = 0 2280 compatibleCoeffs(icoeff2,icoeff) = 0 2281 end if 2282 !!TEST_AM 2283 ! compatibleCoeffs(icoeff,icoeff2) = 0 2284 ! compatibleCoeffs(icoeff2,icoeff) = 0 2285 !!TEST_AM 2286 ! end if 2287 ! end if 2288 2289 ! TODO: hexu: Check if this is only valid for orthogonal structures. As it takes the x component of the rprimd 2290 ! to compare. 2291 ! Here it checks: for two pairs of atom with icoeff (a-b) and icoeff2 (c-d). 2292 !. (a) a-d along x. (b) a-d along y. (c) a-d along z are within the range. 2293 ! As a-b and c-d are within the cutoff, and a==c, there is no more need to check. 2294 2295 2296 if(icoeff<=ncoeff_symsym.and.icoeff2<=ncoeff_symsym)then !Check combination of irreducible bodies and their symmetric equivalent 2297 if(list_symcoeff(2,icoeff,1)/= list_symcoeff(2,icoeff2,1)) then 2298 ABI_BUG("The first components of the pairs in the coefficient are not equivalent.") 2299 end if 2300 2301 if(abs(dist(1,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) & ! rx(a, b) 2302 & -dist(1,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) & ! rx(a, d ) 2303 & >= (rprimd(1,1) + rprimd(1,2) + rprimd(1,3))*sc_size(1) .or. & 2304 abs(dist(2,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) & ! ry(a, b) 2305 & -dist(2,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) & ! ry(a, d) 2306 & >= (rprimd(2,1) + rprimd(2,2) + rprimd(2,3))*sc_size(2) .or. & 2307 abs(dist(3,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff,1),list_symcoeff(4,icoeff,1)) & ! rz(a, b) 2308 & -dist(3,list_symcoeff(2,icoeff,1),list_symcoeff(3,icoeff2,1),list_symcoeff(4,icoeff2,1))) & ! rz(a, d) 2309 & >= (rprimd(3,1) + rprimd(3,2) + rprimd(3,3))*sc_size(3))then 2310 compatibleCoeffs(icoeff,icoeff2) = 0 2311 compatibleCoeffs(icoeff2,icoeff) = 0 2312 endif 2313 endif 2314 end do !end icoeff 2315 end do !icoeff2 2316 2317 ABI_FREE(dist) 2318 ! Compute all the combination of coefficient up to the given order (get the number) 2319 if(need_verbose)then 2320 write(message,'(1a)')' Compute the number of possible combinations' 2321 call wrtout(std_out,message,'COLL') 2322 end if 2323 2324 ABI_MALLOC(list_coeff,(0)) 2325 ABI_MALLOC(list_combination,(0,0)) 2326 icoeff = 1 2327 icoeff2 = 0 2328 nirred_comb = 0 2329 iirred_comb = 0 2330 call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,& 2331 & list_coeff,list_combination,icoeff,max_power_strain,natom,ncoeff_sym,& 2332 & ncoeff_symsym,iirred_comb,nirred_comb,nstr_sym,icoeff,nrpt,nsym,1,power_disps(1),power_disps(2),symbols,comm,& 2333 & nbody=option,compute=.false.,anharmstr=need_anharmstr,spcoupling=need_spcoupling,& 2334 & only_odd_power=need_only_odd_power,only_even_power=need_only_even_power,disp=need_dispterms) 2335 ABI_FREE(list_coeff) 2336 ABI_FREE(list_combination) 2337 2338 ! Output how much we found 2339 if(need_verbose)then 2340 write(message,"(1a,I10)") " -Number of irreducible pairs within cutoff: ", ncoeff_sym 2341 call wrtout(std_out,message,'COLL') 2342 write(message,"(1a,I10)") " -Number of combinations of irreducible pairs: ", nirred_comb 2343 call wrtout(std_out,message,'COLL') 2344 !write(message,"(1a,I10)") " -Number of possible symmetric combinations: ", ncombination 2345 !call wrtout(std_out,message,'COLL') 2346 write(message,'(a,a)') ch10,' Compute the combinations of irreducible pairs' 2347 call wrtout(std_out,message,'COLL') 2348 end if 2349 ABI_MALLOC(list_coeff,(0)) 2350 ABI_MALLOC(list_combination_tmp,(power_disps(2),nirred_comb)) 2351 icoeff = 1 2352 icoeff2 = 0 2353 iirred_comb = 0 2354 list_combination_tmp = 0 2355 ! Compute all the combination of coefficient up to the given order 2356 call computeCombinationFromList(cell,compatibleCoeffs,list_symcoeff,list_symstr,& 2357 & list_coeff,list_combination_tmp,icoeff,max_power_strain,natom,& 2358 & ncoeff_sym,ncoeff_symsym,iirred_comb,nirred_comb,nstr_sym,ncombination,nrpt,nsym,1,power_disps(1),& 2359 & power_disps(2),symbols,comm,nbody=option,compute=.true.,& 2360 & anharmstr=need_anharmstr,spcoupling=need_spcoupling,& 2361 & only_odd_power=need_only_odd_power,only_even_power=need_only_even_power,disp=need_dispterms) 2362 ABI_FREE(list_coeff) 2363 nirred_comb = size(list_combination_tmp,2) 2364 2365 ! If we want to compute equivalent symmetric combinations go here. 2366 if(need_compute_symmetric)then 2367 if(need_verbose)then 2368 write(message,'(1a)')' Distribute irreducible combinations over CPU' 2369 call wrtout(std_out,message,'COLL') 2370 write(message,'(1a)')' Compute symmetric combinations of combinations of irreducible pairs' 2371 call wrtout(std_out,message,'COLL') 2372 write(message,'(3a)')' ---> Try to match number of CPU to number combinations of irreducible pairs',ch10,& 2373 & ' for max. speedup' 2374 call wrtout(std_out,message,'COLL') 2375 endif 2376 2377 ! COUNT IRREDUCIBLE COMBINATIONS FOR EACH PROCESSOR 2378 ncombi_alone = mod(nirred_comb,nproc) 2379 my_ncombi_simple = int(aint(real(nirred_comb,sp)/(nproc))) 2380 if(ncombi_alone == 0 .and. nirred_comb >= nproc)then !ncombi > nproc and no remainder 2381 my_ncombi_start = (my_ncombi_simple * my_rank) + 1 2382 my_ncombi_end = my_ncombi_start + my_ncombi_simple - 1 2383 else if(nirred_comb < nproc)then !ncombi smaller than nproc 2384 if(my_rank + 1 <= nirred_comb)then !myrank smaller than ncombi 2385 my_ncombi_start = my_rank + 1 2386 my_ncombi_end = my_ncombi_start 2387 else 2388 my_ncombi_start = nirred_comb + 1 !myrank bigger than ncombi 2389 my_ncombi_end = nirred_comb + 1 2390 endif 2391 else if(nirred_comb > nproc .and. ncombi_alone /= 0)then !ncombi > nproc and remainder 2392 if(my_rank >= (nproc-ncombi_alone)) then 2393 my_ncombi_start = (my_ncombi_simple * my_rank) + 1 + (my_rank - nproc + ncombi_alone) 2394 my_ncombi_end = my_ncombi_start + my_ncombi_simple 2395 else 2396 my_ncombi_start = (my_ncombi_simple * my_rank) + 1 2397 my_ncombi_end = my_ncombi_start + my_ncombi_simple - 1 2398 endif 2399 end if 2400 2401 if(my_ncombi_end <= nirred_comb)then 2402 my_nirred = my_ncombi_end+1-my_ncombi_start 2403 else 2404 my_nirred = 0 2405 endif 2406 2407 2408 !write(std_out,*) "my_rank", my_rank, "my_ncmobi_start", my_ncombi_start, "my_ncombi_end", my_ncombi_end,'my_nirred',my_nirred 2409 2410 !COPY IRREDUCIBLE COMBINTATIONS TO BE DONE TO EACH PROCESSOR 2411 !my_nirred = my_ncombi_end - my_ncombi_start + 1 2412 ABI_MALLOC(my_list_combination_tmp,(power_disps(2),my_nirred)) 2413 if(my_nirred /= 0)my_list_combination_tmp(:,:) = list_combination_tmp(:,my_ncombi_start:my_ncombi_end) 2414 ABI_MALLOC(my_index_irredcomb,(my_nirred)) 2415 !DEALLOCATE list_combination_tmp 2416 ABI_FREE(list_combination_tmp) 2417 2418 !COUNT SYMMETRIC COMBINATIONS TO IRREDUCIBLE COMBINATIONS ON EACH PROCESSOR 2419 2420 2421 if (my_nirred /= 0 ) then 2422 my_ncombi = my_nirred*(nsym**(power_disps(2))) 2423 end if 2424 2425 !Copy irreducible combinations from list_combination tmp into my_list_combination on rank i 2426 !Make my_list_combination large enough for all symmetric combinations 2427 2428 !ABI_MALLOC(my_list_combination,(power_disps(2),my_ncombi)) 2429 !my_list_combination = 0 2430 2431 2432 2433 !ABI_MALLOC(my_index_irredcomb,(my_ncombi_end-my_ncombi_start+1)) 2434 ! if(my_ncombi /= 0)then 2435 ! do i=1,my_nirred 2436 ! my_list_combination(:,1+(i-1)*(nsym**(power_disps(2)))) = my_list_combination_tmp(:,i) 2437 ! my_index_irredcomb(i) = 1+(i-1)*(nsym**(power_disps(2))) 2438 ! end do 2439 ! endif 2440 2441 !FOREGET ABOUT MY_LIST_COMBINATION_TMP 2442 !COMPUTE SYMMETRIC COMBINATIONS 2443 if(my_ncombi /= 0)then 2444 do i=1,my_nirred 2445 associate(comb => my_list_combination_tmp(:, i)) 2446 !ABI_MALLOC(dummylist,(0)) 2447 ! Get number of strain and displacements for this term 2448 ndisp = 0 2449 nstrain = 0 2450 do ii = 1,power_disps(2) 2451 if(comb(ii) > 0 .and. comb(ii) <= ncoeff_symsym)then 2452 ndisp = ndisp + 1 2453 else if(comb(ii) >= ncoeff_symsym)then 2454 nstrain = nstrain + 1 2455 endif 2456 enddo 2457 compute_sym = .true. 2458 iterm = my_index_irredcomb(i)-1 2459 call computeSymmetricCombinations(my_array_combination,list_symcoeff,list_symstr,ndisp,nsym,& 2460 & comb(:ndisp+nstrain),power_disps(2),& 2461 & my_ncombi,ncoeff_symsym,nstr_sym,nstrain, & 2462 & compatibleCoeffs,compute_sym,comm,only_even=need_only_even_power) 2463 2464 !ABI_FREE(dummylist) 2465 end associate 2466 enddo 2467 endif 2468 2469 ABI_FREE(my_list_combination_tmp) 2470 2471 ! Delete double combinations on each processor 2472 if(need_verbose .and. my_nirred /= 0)then 2473 write(message,'(1a,I4)')' Reduce reducible symmetric combinations on processor: ', my_rank+1 2474 call wrtout(std_out,message,'PERS') 2475 endif 2476 ! TODO: put my_array_combination into my_list_combination. 2477 !call reduce_zero_combinations(my_list_combination) 2478 call my_array_combination%tostatic(my_list_combination, size1=power_disps(2)) 2479 call my_array_combination%finalize() 2480 2481 2482 ! Gather the Results into list_combination_tmp 2483 my_ncombi = size(my_list_combination,2) 2484 ABI_MALLOC(irank_ncombi,(nproc)) 2485 call xmpi_allgather(my_ncombi,irank_ncombi,comm,ierr) 2486 if(need_verbose)then 2487 write(message,'(1a)')' Reduction on all processors finished. Gather results.' 2488 call wrtout(std_out,message,'COLL') 2489 endif 2490 ABI_MALLOC(list_combination_tmp,(power_disps(2),sum(irank_ncombi))) 2491 ABI_MALLOC(offsets,(nproc)) 2492 offsets(1) = 0 2493 do i=1,nproc 2494 offsets(i) = sum(irank_ncombi(:i-1))*power_disps(2) 2495 enddo 2496 2497 list_combination_tmp = 0 2498 ABI_MALLOC(buffsize,(nproc)) 2499 do i = 1,nproc 2500 buffsize(i) = irank_ncombi(i)*power_disps(2) 2501 enddo 2502 call xmpi_gatherv(my_list_combination,size(my_list_combination),list_combination_tmp,buffsize,offsets,master,comm,ierr) 2503 2504 !Deallocation of variables inside need_symmetric 2505 ABI_FREE(buffsize) 2506 ABI_FREE(my_list_combination) 2507 ABI_FREE(my_index_irredcomb) 2508 ABI_FREE(irank_ncombi) 2509 ABI_FREE(offsets) 2510 endif !compute_symmetric 2511 2512 !Deallocation of arrays outside need_symmetric 2513 ABI_FREE(compatibleCoeffs) 2514 2515 if(iam_master)then 2516 call reduce_zero_combinations(list_combination_tmp) 2517 ncombination = size(list_combination_tmp,2) 2518 !ABI_MALLOC(index_irred,(1)) 2519 !index_irred = 1 2520 if(need_spcoupling)then !Check irreducibility of strain-phonon terms 2521 if(need_verbose)then 2522 write(message,'(1a)')' Reduce reducible Strain-Phonon combinations on master' 2523 call wrtout(std_out,message,'COLL') 2524 endif 2525 2526 block 2527 type(IrreducibleCombinations_T) :: irred_combinations 2528 call irred_combinations%init() 2529 do i=1, ncombination 2530 if(any(list_combination_tmp(:,i) > ncoeff_symsym))then 2531 irreducible=irred_combinations%add_irr(list_combination_tmp(:,i), & 2532 & list_symcoeff, list_symstr, ncoeff_symsym, nsym, power_disps(2)) 2533 if(.not. irreducible) list_combination_tmp(:,i) = 0 2534 endif 2535 end do 2536 call reduce_zero_combinations(list_combination_tmp) 2537 ncombination = size(list_combination_tmp,2) 2538 call irred_combinations%free() 2539 end block 2540 endif 2541 2542 2543 ! do i=2,ncombination 2544 ! if(any(list_combination_tmp(:,i) > ncoeff_symsym))then 2545 ! !irreducible = check_irreducibility(list_combination_tmp(:,i),list_combination_tmp(:,:i-1),list_symcoeff,& 2546 ! & ! list_symstr,ncoeff_symsym,nsym,i-1,power_disps(2),index_irred) 2547 ! !if(.not. irreducible) list_combination_tmp(:,i) = 0 2548 ! endif 2549 ! enddo 2550 ! call reduce_zero_combinations(list_combination_tmp) 2551 ! ncombination = size(list_combination_tmp,2) 2552 ! ABI_FREE(index_irred) 2553 end if !iam_master 2554 if(need_verbose)then 2555 write(message,'(1x,I0,1a)') ncombination,' irreducible combinations generated ' 2556 call wrtout(std_out,message,'COLL') 2557 write(message,'(1a)') ' Finished generating irreducible combinations' 2558 call wrtout(std_out,message,'COLL') 2559 endif 2560 2561 !ABI_FREE(xcart) 2562 !ABI_FREE(xred) 2563 2564 !MPI 2565 if(need_verbose .and. nproc > 1)then 2566 write(message,'(a,a)') ch10,' Redistribute the combinations over the CPU' 2567 call wrtout(std_out,message,'COLL') 2568 end if 2569 2570 call xmpi_bcast(ncombination, master, comm, ierr) 2571 2572 ncoeff_alone = mod(ncombination,nproc) 2573 my_ncoeff = int(aint(real(ncombination,sp)/(nproc))) 2574 2575 if(my_rank >= (nproc-ncoeff_alone)) then 2576 my_ncoeff = my_ncoeff + 1 2577 end if 2578 2579 !Set the buffsize for mpi scatterv 2580 ABI_MALLOC(buffsize,(nproc)) 2581 ABI_MALLOC(buffdispl,(nproc)) 2582 do ii = 1,nproc 2583 buffsize(ii) = int(aint(real(ncombination,sp)/(nproc))*power_disps(2)) 2584 if(ii > (nproc-ncoeff_alone)) then 2585 buffsize(ii) = buffsize(ii) + power_disps(2) 2586 end if 2587 end do 2588 2589 buffdispl(1) = 0 2590 do ii = 2,nproc 2591 buffdispl(ii) = buffdispl(ii-1) + buffsize(ii-1) 2592 end do 2593 2594 ABI_MALLOC(list_combination,(power_disps(2),my_ncoeff)) 2595 list_combination = 0 2596 2597 my_size = my_ncoeff*power_disps(2) 2598 call xmpi_scatterv(list_combination_tmp,buffsize,buffdispl,list_combination,my_size,master,& 2599 & comm,ierr) 2600 2601 ABI_FREE(buffdispl) 2602 ABI_FREE(buffsize) 2603 ABI_FREE(list_combination_tmp) 2604 2605 !write(std_out,*) "DEBUG shape list_symcoeff(:,:,:) before termsfromlist:", shape(list_symcoeff) 2606 !Compute the coefficients from the list of combination 2607 if(need_verbose .and. nproc > 1)then 2608 write(message,'(1a)')' Compute the coefficients' 2609 call wrtout(std_out,message,'COLL') 2610 end if 2611 ABI_MALLOC(coeffs_tmp,(my_ncoeff)) 2612 nterm = nsym 2613 ndisp_max = power_disps(2) 2614 ncoeff_max = my_ncoeff 2615 do ii=1,my_ncoeff 2616 ABI_MALLOC(terms,(nterm)) 2617 !write(std_out,*) "DEBUG list_combination(:,",ii,"): ", list_combination(:,ii) 2618 block 2619 logical :: reverse(ndisp_max) 2620 reverse=.False. 2621 call generateTermsFromList(cell,list_combination(:,ii),list_symcoeff,list_symstr,ncoeff_symsym,& 2622 & ndisp_max,nrpt,nstr_sym,nsym,nterm,terms, reverse=reverse) 2623 call polynomial_coeff_init(one,nterm,coeffs_tmp(ii),terms(1:nterm),check=.true.) 2624 end block 2625 ! Free the terms array 2626 do iterm=1,nterm 2627 call polynomial_term_free(terms(iterm)) 2628 end do 2629 ABI_FREE(terms) 2630 end do 2631 ABI_FREE(cell) 2632 2633 2634 ABI_FREE(list_combination) 2635 ABI_FREE(list_symcoeff) 2636 ABI_FREE(list_symstr) 2637 2638 !Final tranfert 2639 !1- Count the total number of coefficient 2640 ncoeff = 0 2641 do icoeff=1,ncoeff_max 2642 if (abs(coeffs_tmp(icoeff)%coefficient) >tol16) then 2643 ncoeff = ncoeff + 1 2644 end if 2645 end do 2646 2647 !Get the total number of coefficients 2648 !ncoeff_max is the number of total coefficients before the symetries check 2649 !ncoeff_tot is the number of total coefficients after the symetries check 2650 ncoeff_tot = ncoeff!set the output 2651 call xmpi_sum(ncoeff_tot,comm,ierr) 2652 call xmpi_sum(ncoeff_max,comm,ierr) 2653 2654 2655 !Need to redistribute the coefficients over the CPU 2656 !Get the list with the number of coeff on each CPU 2657 !In order to be abble to compute the my_coeffindexes array which is for example: 2658 ! if CPU0 has 200 Coeff and CPU1 has 203 Coeff then 2659 ! for CPU0:my_coeffindexes=>1-200 and for CPU1:my_coeffindexes=>201-403 2660 if(need_verbose .and. nproc > 1)then 2661 write(message,'(1a)')' Redistribute the coefficients over the CPU' 2662 call wrtout(std_out,message,'COLL') 2663 end if 2664 2665 ABI_MALLOC(buffdispl,(nproc)) 2666 buffdispl = 0 2667 buffdispl(my_rank+1) = my_ncoeff 2668 call xmpi_sum(buffdispl,comm,ierr) 2669 ABI_MALLOC(my_coeffindexes,(my_ncoeff)) 2670 ABI_MALLOC(my_coefflist,(my_ncoeff)) 2671 my_coeffindexes = 0 2672 my_coefflist = 0 2673 do icoeff=1,my_ncoeff 2674 my_coefflist(icoeff) = icoeff 2675 if(my_rank==0) then 2676 my_coeffindexes(icoeff) = icoeff 2677 else 2678 my_coeffindexes(icoeff) = sum(buffdispl(1:my_rank)) + icoeff 2679 end if 2680 end do 2681 ABI_FREE(buffdispl) 2682 2683 !Compute the new number of coefficient per CPU 2684 if(need_distributed) then 2685 ncoeff_alone = mod(ncoeff_tot,nproc) 2686 my_newncoeff = int(aint(real(ncoeff_tot,sp)/(nproc))) 2687 if(my_rank >= (nproc-ncoeff_alone)) then 2688 my_newncoeff = my_newncoeff + 1 2689 end if 2690 else 2691 my_newncoeff = ncoeff_tot 2692 end if 2693 2694 ncoeff = my_newncoeff ! Set the output 2695 2696 !2:compute the number of coefficients and the list of the corresponding 2697 ! coefficients for each CPU. 2698 ABI_MALLOC(my_newcoeffindexes,(my_newncoeff)) 2699 if(need_distributed) then 2700 do icoeff=1,my_newncoeff 2701 if(my_rank >= (nproc-ncoeff_alone))then 2702 my_newcoeffindexes(icoeff)=int(aint(real(ncoeff_tot,sp)/(nproc)))*(my_rank)+& 2703 & (my_rank - (nproc-ncoeff_alone)) + icoeff 2704 else 2705 my_newcoeffindexes(icoeff)=(my_newncoeff)*(my_rank) + icoeff 2706 end if 2707 end do 2708 else 2709 do icoeff=1,my_newncoeff 2710 my_newcoeffindexes(icoeff) = icoeff 2711 end do 2712 end if 2713 2714 !2- Transfer 2715 if(.not.need_distributed)then 2716 if(.not.allocated(coefficients))then 2717 ABI_MALLOC(coefficients,(my_newncoeff)) 2718 end if 2719 end if 2720 icoeff = 0! icoeff is the current index in the total list of coefficients 2721 icoeff2 = 0! icoeff2 is the current index in the output coefficients array on each CPU 2722 icoeff3 = 0! icoeff3 is the current index in total new list of coefficients 2723 rank_to_send_save = 0 2724 2725 do icoeff=1,ncoeff_max 2726 ! Need to send the rank with the chosen coefficient 2727 rank_to_send = 0 2728 my_icoeff = 0 2729 do ii=1,my_ncoeff 2730 if (my_coeffindexes(ii)==icoeff) then 2731 my_icoeff = ii 2732 if (abs(coeffs_tmp(my_icoeff)%coefficient) > tol16)then 2733 rank_to_send = my_rank 2734 else 2735 rank_to_send = -1 2736 ! Free the coefficient 2737 call polynomial_coeff_free(coeffs_tmp(ii)) 2738 end if 2739 exit 2740 end if 2741 end do 2742 call xmpi_sum(rank_to_send, comm, ierr) 2743 ! This coefficient is not compute 2744 if (rank_to_send == -1) cycle 2745 2746 ! increase icoeff3 2747 icoeff3 = icoeff3 + 1 2748 2749 ! Find the receiver CPU 2750 rank_to_receive = 0 2751 do ii=1,my_newncoeff 2752 if (my_newcoeffindexes(ii)==icoeff3) then 2753 rank_to_receive = my_rank 2754 end if 2755 end do 2756 call xmpi_sum(rank_to_receive, comm, ierr) 2757 2758 if(need_distributed.and.rank_to_send /= rank_to_send_save) then 2759 if(my_rank == rank_to_send_save)then 2760 ABI_FREE(coeffs_tmp)!Free memory if the current CPU has already distribute 2761 !all its own coefficients 2762 end if 2763 rank_to_send_save = rank_to_send 2764 end if 2765 2766 if(need_distributed.and.my_rank == rank_to_receive)then 2767 if(.not.allocated(coefficients))then 2768 ABI_MALLOC(coefficients,(my_newncoeff)) 2769 end if 2770 end if 2771 2772 2773 if (need_distributed)then 2774 if(my_rank==rank_to_send)then 2775 if(any(my_newcoeffindexes(:)==icoeff3))then 2776 icoeff2 = icoeff2 + 1 2777 ! Get the name of this coefficient 2778 call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.) 2779 call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),& 2780 & coeffs_tmp(my_icoeff)%terms,name=name,check=.false.) 2781 else 2782 call polynomial_coeff_MPIsend(coeffs_tmp(my_icoeff), icoeff, rank_to_receive, comm) 2783 end if 2784 ! Free the coefficient 2785 call polynomial_coeff_free(coeffs_tmp(my_icoeff)) 2786 else 2787 if(any(my_newcoeffindexes(:)==icoeff3))then 2788 icoeff2 = icoeff2 + 1 2789 call polynomial_coeff_MPIrecv(coefficients(icoeff2), icoeff, rank_to_send, comm) 2790 call polynomial_coeff_getName(name,coefficients(icoeff2),symbols,recompute=.TRUE.) 2791 call polynomial_coeff_SetName(name,coefficients(icoeff2)) 2792 end if 2793 end if 2794 else 2795 icoeff2 = icoeff2 + 1 2796 ! Get the name of this coefficient 2797 if(my_rank==rank_to_send)then 2798 call polynomial_coeff_getName(name,coeffs_tmp(my_icoeff),symbols,recompute=.TRUE.) 2799 call polynomial_coeff_init(one,coeffs_tmp(my_icoeff)%nterm,coefficients(icoeff2),& 2800 & coeffs_tmp(my_icoeff)%terms,name=name,check=.false.) 2801 ! Free the coefficient 2802 call polynomial_coeff_free(coeffs_tmp(my_icoeff)) 2803 end if 2804 call polynomial_coeff_broadcast(coefficients(icoeff2),rank_to_send, comm) 2805 end if 2806 end do 2807 2808 !Debug write xml 2809 ! write (filename, "(A9,I2,A4)") "terms_set", my_rank+1,".xml" 2810 ! call polynomial_coeff_writeXML(coefficients,my_newncoeff,filename=filename) 2811 2812 if(need_verbose)then 2813 write(message,'(1x,I0,2a)') ncoeff_tot,' coefficients generated ',ch10 2814 call wrtout(ab_out,message,'COLL') 2815 call wrtout(std_out,message,'COLL') 2816 end if 2817 2818 2819 !Final deallocation 2820 ABI_FREE(symbols) 2821 ABI_FREE(my_coeffindexes) 2822 ABI_FREE(my_newcoeffindexes) 2823 ABI_FREE(my_coefflist) 2824 ABI_SFREE(coeffs_tmp) 2825 2826 end subroutine polynomial_coeff_getNorder
m_polynomial_coeff/polynomial_coeff_getOrder1 [ Functions ]
NAME
polynomial_coeff_getOrder1
FUNCTION
Compute the first order polynomial coefficients from the list
INPUTS
cell(3,nrpt) = indexes of the cells into the supercell (-1 -1 -1, 0 0 0 ...) cutoff_in = cut-off for the inter atomic forces constants list_symcoeff(6,ncoeff_sym,nsym) = array with the list of the coefficients, for each coefficients (ncoeff_sym), we store the symmetrics(nsym) the 6th first dimensions are : 1 = direction of the IFC 2 = index of the atom number 1 (1=>natom) 3 = index of the atom number 2 (1=>natom) 4 = indexes of the cell of the second atom (the atom number 1 is always in the cell 0 0 0) 5 = weight of the term (-1 or 1) 6 = indexes of the symmetric natom = number of atoms in the unit cell nrpt = number of cell nsym = number of symmetries in the system symbols(natom) = array with the symbols of each atoms (Sr,O,Ti,...) comm = MPI communicator
OUTPUT
polynomial_coeff<(type(polynomial_coeff_type)>(ncoeff_out) = array of datatype with the polynomial_coeff ncoeff_out = number of coefficients
SOURCE
3754 subroutine polynomial_coeff_getOrder1(cell,coeffs_out,list_symcoeff,& 3755 & natom,ncoeff_out,ncoeff,nrpt,nsym,& 3756 & symbols) 3757 3758 implicit none 3759 3760 !Arguments ------------------------------------ 3761 !scalars 3762 integer,intent(in) :: natom,ncoeff,nsym,nrpt 3763 integer,intent(out) :: ncoeff_out 3764 !arrays 3765 integer,intent(in) :: cell(3,nrpt) 3766 integer,intent(in) :: list_symcoeff(6,ncoeff,nsym) 3767 character(len=5),intent(in) :: symbols(natom) 3768 type(polynomial_coeff_type),allocatable,intent(inout) :: coeffs_out(:) 3769 !Local variables------------------------------- 3770 !scalar 3771 integer :: ia,ib,icoeff,icoeff_tmp,irpt,irpt_ref 3772 integer :: isym,iterm,mu,ncoeff_max,ndisp,nstrain,nterm_max 3773 real(dp):: coefficient,weight 3774 !arrays 3775 integer,allocatable :: atindx(:,:),cells(:,:,:),dir_int(:) 3776 integer,allocatable :: power_disps(:),power_strain(:),strain(:) 3777 character(len=1) :: mutodir(9) = (/"x","y","z","1","2","3","4","5","6"/) 3778 character(len=200):: name 3779 character(len=500) :: message 3780 type(polynomial_term_type),dimension(:),allocatable :: terms 3781 type(polynomial_coeff_type),allocatable :: coeffs_tmp(:) 3782 !TEST_AM 3783 character(len=fnlen) :: filename 3784 !TEST_AM 3785 ! ************************************************************************* 3786 3787 !Initialisation of variables 3788 nterm_max = nsym 3789 ncoeff_max = ncoeff 3790 ndisp = 1 3791 nstrain = 0 3792 ABI_MALLOC(coeffs_tmp,(ncoeff_max)) 3793 ABI_MALLOC(terms,(nterm_max)) 3794 3795 3796 icoeff_tmp = 0 3797 ABI_MALLOC(atindx,(2,ndisp)) 3798 ABI_MALLOC(cells,(3,2,ndisp)) 3799 ABI_MALLOC(dir_int,(ndisp)) 3800 ABI_MALLOC(power_disps,(ndisp)) 3801 ABI_MALLOC(power_strain,(nstrain)) 3802 ABI_MALLOC(strain,(nstrain)) 3803 3804 !Found the ref cell 3805 irpt_ref = 1 3806 do irpt=1,nrpt 3807 if(all(cell(:,irpt)==0))then 3808 irpt_ref = irpt 3809 exit 3810 end if 3811 end do 3812 3813 write(message,'(3a)') " Irreductible coefficient and associated atom 1, atom 2 and direction:",ch10,& 3814 & " for the 1st order" 3815 call wrtout(std_out,message,'COLL') 3816 3817 do icoeff=1,ncoeff 3818 ! Reset counter 3819 iterm = 0 3820 coefficient = one 3821 do isym=1,nsym 3822 ndisp = 1 3823 nstrain = 0 3824 mu = list_symcoeff(1,icoeff,isym) 3825 ia = list_symcoeff(2,icoeff,isym) 3826 ib = list_symcoeff(3,icoeff,isym) 3827 irpt = list_symcoeff(4,icoeff,isym) 3828 weight = list_symcoeff(5,icoeff,isym) 3829 ! Fill First term arrays 3830 atindx(1,1) = ia; atindx(2,1) = ib; 3831 dir_int(1) = mu 3832 power_disps(1) = 1 3833 cells(:,1,1) = (/0,0,0/) 3834 cells(:,2,1) = cell(:,irpt) 3835 iterm = iterm + 1 3836 call polynomial_term_init(atindx,cells,dir_int,ndisp,nstrain,terms(iterm),& 3837 & power_disps,power_strain,strain,weight,check=.true.) 3838 end do!end do sym 3839 3840 if(iterm > 0)then 3841 ! increase coefficients and set it 3842 icoeff_tmp = icoeff_tmp + 1 3843 call polynomial_coeff_init(coefficient,iterm,coeffs_tmp(icoeff_tmp),terms(1:iterm),check=.true.) 3844 end if 3845 3846 ! Deallocate the terms 3847 do iterm=1,nterm_max 3848 call polynomial_term_free(terms(iterm)) 3849 end do 3850 end do!end do coeff_sym 3851 3852 ABI_FREE(terms) 3853 ABI_FREE(atindx) 3854 ABI_FREE(cells) 3855 ABI_FREE(dir_int) 3856 ABI_FREE(power_disps) 3857 ABI_FREE(power_strain) 3858 ABI_FREE(strain) 3859 3860 !Count the number of terms 3861 ncoeff_out = 0 3862 do icoeff_tmp=1,ncoeff_max 3863 if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then 3864 ncoeff_out = ncoeff_out + 1 3865 end if 3866 end do 3867 3868 !Transfer in the final array 3869 ABI_MALLOC(coeffs_out,(ncoeff_out)) 3870 icoeff = 0 3871 do icoeff_tmp=1,ncoeff_max 3872 if (abs(coeffs_tmp(icoeff_tmp)%coefficient) > tol16)then 3873 ! Get the name of this coefficient 3874 call polynomial_coeff_getName(name,coeffs_tmp(icoeff_tmp),symbols,recompute=.TRUE.) 3875 ! Increase icoeff and fill the coeffs_out array 3876 icoeff = icoeff + 1 3877 call polynomial_coeff_init(one,coeffs_tmp(icoeff_tmp)%nterm,& 3878 & coeffs_out(icoeff),coeffs_tmp(icoeff_tmp)%terms,name=name) 3879 3880 write(message,'(2a)')' ',trim(name) 3881 call wrtout(std_out,message,'COLL') 3882 3883 do iterm = 1,coeffs_tmp(icoeff_tmp)%nterm 3884 write(message,'(a,I0,a,I0,2a)') ' Atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(1,1),& 3885 & ' and atom ',coeffs_tmp(icoeff_tmp)%terms(iterm)%atindx(2,1),& 3886 & ' in the direction ',mutodir(coeffs_tmp(icoeff_tmp)%terms(iterm)%direction(1)) 3887 if(any(coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(:,2,1)/=0))then 3888 write(message,'(2a,I0,a,I0,a,I0,a)') trim(message),' in the cell ',& 3889 & coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(1,2,1),' ',& 3890 & coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(2,2,1),' ',& 3891 & coeffs_tmp(icoeff_tmp)%terms(iterm)%cell(3,2,1),'.' 3892 end if 3893 call wrtout(std_out,message,'COLL') 3894 end do 3895 end if 3896 end do 3897 3898 !TEST_AM 3899 filename = "terms_1st_order.xml" 3900 call polynomial_coeff_writeXML(coeffs_out,ncoeff_out,filename=filename) 3901 !TEST_AM 3902 3903 write(message,'(a,1x,I0,a)') ch10,& 3904 & ncoeff_out,' fitted coefficients for the 1st order ' 3905 call wrtout(ab_out,message,'COLL') 3906 call wrtout(std_out,message,'COLL') 3907 3908 !Deallocation 3909 do icoeff=1,ncoeff_max 3910 call polynomial_coeff_free(coeffs_tmp(icoeff)) 3911 end do 3912 ABI_FREE(coeffs_tmp) 3913 3914 end subroutine polynomial_coeff_getOrder1
m_polynomial_coeff/polynomial_coeff_init [ Functions ]
NAME
polynomial_coeff_init
FUNCTION
Initialize a polynomial_coeff datatype
INPUTS
name = Name of the polynomial_coeff (Sr_y-O1_y)^3) for example nterm = Number of terms (short range interaction) for this polynomial_coeff coefficient = Value of the coefficient of this term terms(nterm)<type(polynomial_term)> = array of polynomial_term_type check = TRUE if this list of terms has to be check. We remove the symetric of equivalent terms for example: ((Sr_y-O1_y)^1 and -1*(Sr_y-O1_y)^1 => zero FALSE, defaut, do nothing
OUTPUT
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype to be initialized
SOURCE
196 subroutine polynomial_coeff_init(coefficient,nterm,polynomial_coeff,terms,name,check) 197 198 implicit none 199 200 !Arguments ------------------------------------ 201 !scalars 202 integer, intent(in) :: nterm 203 real(dp),intent(in) :: coefficient 204 logical,optional,intent(in) :: check 205 !arrays 206 character(len=200),optional,intent(in) :: name 207 type(polynomial_term_type),intent(in) :: terms(nterm) 208 type(polynomial_coeff_type), intent(out) :: polynomial_coeff 209 !Local variables------------------------------- 210 !scalar 211 integer :: iterm1,iterm2 212 integer :: ii,nterm_tmp 213 real(dp):: coefficient_tmp 214 logical :: check_in 215 !arrays 216 real(dp) :: weights(nterm) 217 character(len=200) :: name_tmp 218 ! ************************************************************************* 219 !First free before initilisation 220 call polynomial_coeff_free(polynomial_coeff) 221 check_in = .false. 222 if(present(check)) check_in = check 223 if(check_in)then 224 ! Check if the list of term is available or contains identical terms 225 ! in this case, remove all the not needed terms 226 nterm_tmp = 0 227 weights(:) = one 228 do iterm1=1,nterm 229 if(abs(weights(iterm1)) < tol16)cycle ! FIXME: do nothing? 230 weights(iterm1) = terms(iterm1)%weight 231 do iterm2=iterm1+1,nterm 232 if(abs(weights(iterm2)) < tol16)cycle 233 ! if the terms are identical we check the weight 234 if(terms(iterm1)==terms(iterm2))then 235 weights(iterm1) = weights(iterm1) + terms(iterm2)%weight 236 weights(iterm2) = 0 237 end if 238 end do 239 if(abs(weights(iterm1)) > tol16) then 240 ! FIXME: weights(iterm1) = 1 ? 241 weights(iterm1)= anint(weights(iterm1)/weights(iterm1)) 242 end if 243 end do 244 245 246 ! Count the number of terms 247 ! TODO: replace by: ?? 248 !nterm_tmp=count(abs(weights) > tol16) 249 nterm_tmp = 0 250 do iterm1=1,nterm 251 if(abs(weights(iterm1)) > tol16)then 252 nterm_tmp = nterm_tmp + 1 253 end if 254 end do 255 256 if (nterm_tmp ==0)then 257 coefficient_tmp = 0.0 258 else 259 coefficient_tmp = coefficient 260 end if 261 262 else 263 nterm_tmp = nterm 264 coefficient_tmp = coefficient 265 weights(:) = terms(:)%weight 266 end if!end Check 267 268 if(present(name))then 269 name_tmp = name 270 else 271 name_tmp = "" 272 end if 273 274 275 !Initilisation 276 polynomial_coeff%name = name_tmp 277 polynomial_coeff%nterm = nterm_tmp 278 polynomial_coeff%coefficient = coefficient_tmp 279 ABI_MALLOC(polynomial_coeff%terms,(polynomial_coeff%nterm)) 280 iterm1 = 0 281 do ii = 1,nterm 282 if(abs(weights(ii)) > tol16)then 283 iterm1 = iterm1 + 1 284 call polynomial_term_init(terms(ii)%atindx,terms(ii)%cell,terms(ii)%direction,terms(ii)%ndisp,& 285 & terms(ii)%nstrain,polynomial_coeff%terms(iterm1),terms(ii)%power_disp,& 286 & terms(ii)%power_strain,terms(ii)%strain,terms(ii)%weight) 287 end if 288 end do 289 end subroutine polynomial_coeff_init
m_polynomial_coeff/polynomial_coeff_list_free [ Functions ]
NAME
polynomial_coeff_list_free
FUNCTION
Free polynomial_coeff datatype
INPUTS
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
OUTPUT
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
SOURCE
353 subroutine polynomial_coeff_list_free(polynomial_coeff_list) 354 355 implicit none 356 357 !Arguments ------------------------------------ 358 !scalars 359 !arrays 360 type(polynomial_coeff_type),allocatable, intent(inout) :: polynomial_coeff_list(:) 361 !Local variables------------------------------- 362 !scalar 363 integer :: i,ncoeff 364 !arrays 365 366 ! ************************************************************************* 367 368 !Free output 369 if(allocated(polynomial_coeff_list))then 370 ncoeff = size(polynomial_coeff_list) 371 do i=1,ncoeff 372 call polynomial_coeff_free(polynomial_coeff_list(i)) 373 enddo 374 ABI_FREE(polynomial_coeff_list) 375 endif 376 377 end subroutine polynomial_coeff_list_free
m_polynomial_coeff/polynomial_coeff_MPIrecv [ Functions ]
NAME
polynomial_coeff_MPIrecv
FUNCTION
MPI receive the polynomial_coefficent datatype
INPUTS
tag = tag of the message to receive source = rank of Source comm = MPI communicator
SIDE EFFECTS
coefficients<type(polynomial_coefficent_type)>= polynomial_coeff datatype
SOURCE
755 subroutine polynomial_coeff_MPIrecv(coefficients, tag, source, comm) 756 757 implicit none 758 759 !Arguments ------------------------------------ 760 !array 761 type(polynomial_coeff_type),intent(inout) :: coefficients 762 integer, intent(in) :: source,comm,tag 763 764 !Local variables------------------------------- 765 !scalars 766 integer :: ierr,ii 767 !arrays 768 769 ! ************************************************************************* 770 771 if (xmpi_comm_size(comm) == 1) return 772 773 774 ! Free the output 775 call polynomial_coeff_free(coefficients) 776 777 ! Transmit variables 778 call xmpi_recv(coefficients%name, source, 9*tag+0, comm, ierr) 779 call xmpi_recv(coefficients%nterm, source, 9*tag+1, comm, ierr) 780 call xmpi_recv(coefficients%coefficient, source, 9*tag+2, comm, ierr) 781 782 !Allocate arrays on the other nodes. 783 ABI_MALLOC(coefficients%terms,(coefficients%nterm)) 784 do ii=1,coefficients%nterm 785 call polynomial_term_free(coefficients%terms(ii)) 786 end do 787 788 ! Set the number of term on each node (needed for allocations of array) 789 do ii = 1,coefficients%nterm 790 call xmpi_recv(coefficients%terms(ii)%ndisp, source, 9*tag+3, comm, ierr) 791 call xmpi_recv(coefficients%terms(ii)%nstrain, source, 9*tag+4, comm, ierr) 792 end do 793 794 ! Allocate arrays on the other nodes 795 do ii = 1,coefficients%nterm 796 ABI_MALLOC(coefficients%terms(ii)%atindx,(2,coefficients%terms(ii)%ndisp)) 797 coefficients%terms(ii)%atindx = 0 798 ABI_MALLOC(coefficients%terms(ii)%direction,(coefficients%terms(ii)%ndisp)) 799 ABI_MALLOC(coefficients%terms(ii)%cell,(3,2,coefficients%terms(ii)%ndisp)) 800 ABI_MALLOC(coefficients%terms(ii)%power_disp,(coefficients%terms(ii)%ndisp)) 801 ABI_MALLOC(coefficients%terms(ii)%power_strain,(coefficients%terms(ii)%nstrain)) 802 ABI_MALLOC(coefficients%terms(ii)%strain,(coefficients%terms(ii)%nstrain)) 803 end do 804 805 ! Transfert value 806 do ii = 1,coefficients%nterm 807 call xmpi_recv(coefficients%terms(ii)%weight, source, 9*tag+5, comm, ierr) 808 call xmpi_recv(coefficients%terms(ii)%atindx, source, 9*tag+6, comm, ierr) 809 call xmpi_recv(coefficients%terms(ii)%direction, source, 9*tag+7, comm, ierr) 810 call xmpi_recv(coefficients%terms(ii)%cell, source, 9*tag+8, comm, ierr) 811 call xmpi_recv(coefficients%terms(ii)%power_disp, source, 9*tag+9, comm, ierr) 812 call xmpi_recv(coefficients%terms(ii)%power_strain, source, 9*tag+10, comm, ierr) 813 call xmpi_recv(coefficients%terms(ii)%strain, source, 9*tag+11, comm, ierr) 814 end do 815 816 end subroutine polynomial_coeff_MPIrecv
m_polynomial_coeff/polynomial_coeff_MPIsend [ Functions ]
NAME
polynomial_coeff_MPIsend
FUNCTION
MPI send the polynomial_coefficent datatype
INPUTS
tag = tag of the message to send dest= rank of Dest comm= MPI communicator
SIDE EFFECTS
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
SOURCE
693 subroutine polynomial_coeff_MPIsend(coefficients, tag, dest, comm) 694 695 implicit none 696 697 !Arguments ------------------------------------ 698 !array 699 type(polynomial_coeff_type),intent(inout) :: coefficients 700 integer, intent(in) :: dest,comm,tag 701 702 !Local variables------------------------------- 703 !scalars 704 integer :: ierr,ii 705 integer :: my_rank 706 !arrays 707 708 ! ************************************************************************* 709 710 if (xmpi_comm_size(comm) == 1) return 711 712 my_rank = xmpi_comm_rank(comm) 713 ! Transmit variables 714 call xmpi_send(coefficients%name, dest, 9*tag+0, comm, ierr) 715 call xmpi_send(coefficients%nterm, dest, 9*tag+1, comm, ierr) 716 call xmpi_send(coefficients%coefficient, dest, 9*tag+2, comm, ierr) 717 718 ! Set the number of term on each node (needed for allocations of array) 719 do ii = 1,coefficients%nterm 720 call xmpi_send(coefficients%terms(ii)%ndisp, dest, 9*tag+3, comm, ierr) 721 call xmpi_send(coefficients%terms(ii)%nstrain, dest, 9*tag+4, comm, ierr) 722 end do 723 724 ! Transfert value 725 do ii = 1,coefficients%nterm 726 call xmpi_send(coefficients%terms(ii)%weight, dest, 9*tag+5, comm, ierr) 727 call xmpi_send(coefficients%terms(ii)%atindx, dest, 9*tag+6, comm, ierr) 728 call xmpi_send(coefficients%terms(ii)%direction, dest, 9*tag+7, comm, ierr) 729 call xmpi_send(coefficients%terms(ii)%cell, dest, 9*tag+8, comm, ierr) 730 call xmpi_send(coefficients%terms(ii)%power_disp, dest, 9*tag+9, comm, ierr) 731 call xmpi_send(coefficients%terms(ii)%power_strain, dest, 9*tag+10, comm, ierr) 732 call xmpi_send(coefficients%terms(ii)%strain, dest, 9*tag+11, comm, ierr) 733 end do 734 735 end subroutine polynomial_coeff_MPIsend
m_polynomial_coeff/polynomial_coeff_setCoefficient [ Functions ]
NAME
polynomial_coeff_setCoefficient
FUNCTION
set the coefficient for of polynomial_coeff
INPUTS
coefficient = coefficient of this coefficient
OUTPUT
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
SOURCE
398 subroutine polynomial_coeff_setCoefficient(coefficient,polynomial_coeff) 399 400 implicit none 401 402 !Arguments ------------------------------------ 403 !scalars 404 real(dp),intent(in) :: coefficient 405 !arrays 406 type(polynomial_coeff_type), intent(inout) :: polynomial_coeff 407 !Local variables------------------------------- 408 !scalar 409 !arrays 410 ! ************************************************************************* 411 412 polynomial_coeff%coefficient = coefficient 413 414 end subroutine polynomial_coeff_setCoefficient
m_polynomial_coeff/polynomial_coeff_setName [ Functions ]
NAME
polynomial_coeff_setName
FUNCTION
set the name of a polynomial_coeff type
INPUTS
name = name of the coeff
OUTPUT
polynomial_coeff<type(polynomial_coeff)> = polynomial_coeff datatype
SOURCE
433 subroutine polynomial_coeff_setName(name,polynomial_coeff) 434 435 implicit none 436 437 !Arguments ------------------------------------ 438 !scalars 439 !arrays 440 character(len=200),intent(in) :: name 441 type(polynomial_coeff_type), intent(inout) :: polynomial_coeff 442 !Local variables------------------------------- 443 !scalar 444 !arrays 445 ! ************************************************************************* 446 447 polynomial_coeff%name = name 448 449 end subroutine polynomial_coeff_setName
m_polynomial_coeff/polynomial_coeff_type [ Types ]
NAME
polynomial_coeff_type
FUNCTION
structure for a polynomial coefficient contains the value of the coefficient and a list of terms (displacements and/or strain) relating to the coefficient by symmetry
SOURCE
133 type, public :: polynomial_coeff_type 134 135 character(len=200) :: name = "" 136 ! Name of the polynomial_coeff (Sr_y-O1_y)^3) for example 137 138 integer :: nterm = 0 139 ! Number of terms (short range interaction) for this polynomial_coeff 140 141 real(dp) :: coefficient = zero 142 ! coefficient = value of the coefficient of this term 143 ! \frac{\partial E^{k}}{\partial \tau^{k}} 144 145 type(polynomial_term_type),dimension(:),allocatable :: terms 146 ! polynomial_term(nterm)<type(polynomial_term)> 147 ! contains all the displacements for this coefficient 148 149 end type polynomial_coeff_type
m_polynomial_coeff/reduce_zero_combinations [ Functions ]
NAME
reduce_zero_combinations
FUNCTION
Sort out list of zeros in a list of integers
INPUTS
OUTPUT
SOURCE
4472 subroutine reduce_zero_combinations(combination_list) 4473 !Arguments ------------------------------------ 4474 implicit none 4475 4476 !Arguments ------------------------------------ 4477 !scalar 4478 !array 4479 integer,allocatable,intent(inout) :: combination_list(:,:) 4480 !local 4481 !variable 4482 integer :: i,j 4483 integer,allocatable :: combination_list_tmp(:,:) 4484 !array 4485 ! ************************************************************************* 4486 4487 !Reduce zero combinations 4488 i = 0 4489 do j=1,size(combination_list,2) 4490 if(any(combination_list(:,j) /= 0))then 4491 i = i + 1 4492 endif 4493 enddo 4494 4495 ABI_MALLOC(combination_list_tmp,(size(combination_list,1),i)) 4496 i = 0 4497 do j=1,size(combination_list,2) 4498 if(any(combination_list(:,j) /= 0))then 4499 i = i + 1 4500 combination_list_tmp(:,i) = combination_list(:,j) 4501 endif 4502 enddo 4503 ABI_FREE(combination_list) 4504 ABI_MALLOC(combination_list,(size(combination_list_tmp,1),i)) 4505 combination_list = combination_list_tmp 4506 ABI_FREE(combination_list_tmp) 4507 4508 end subroutine reduce_zero_combinations
m_polynomial_coeff/sort_combination [ Functions ]
NAME
sort_combination
FUNCTION
Sort a list of integer from small to large if it contains zeros will be put to highest indexe
INPUTS
OUTPUT
SOURCE
4260 subroutine sort_combination(combination,n_int) 4261 !Arguments ------------------------------------ 4262 implicit none 4263 4264 !Arguments ------------------------------------ 4265 !scalar 4266 integer,intent(in) :: n_int 4267 !array 4268 integer,intent(inout) :: combination(n_int) 4269 !local 4270 !variable 4271 integer :: j,k,cnt,tmp_int1,tmp_int2 4272 !array 4273 ! ************************************************************************* 4274 4275 j=2 4276 do while(j <= n_int) 4277 k = j 4278 cnt = 1 4279 do while(k >= 2 .and. cnt == 1) 4280 if(combination(k-1) > combination(k) .and. combination(k) > 0)then 4281 tmp_int1 = combination(k-1) 4282 tmp_int2 = combination(k) 4283 combination(k) = tmp_int1 4284 combination(k-1) = tmp_int2 4285 k=k-1 4286 else 4287 cnt = cnt + 1 4288 end if 4289 end do 4290 j = j+1 4291 end do 4292 4293 end subroutine sort_combination
m_polynomial_coeff/sort_combination_list [ Functions ]
NAME
sort_combination_list
FUNCTION
Sort a list of integer list from small to large if it contains zeros will be put to highest index
INPUTS
OUTPUT
SOURCE
4310 subroutine sort_combination_list(combination_list,n_int,n_list) 4311 !Arguments ------------------------------------ 4312 implicit none 4313 4314 !Arguments ------------------------------------ 4315 !scalar 4316 integer,intent(in) :: n_int,n_list 4317 !array 4318 integer,intent(inout) :: combination_list(n_int,n_list) 4319 !local 4320 !variable 4321 integer :: i 4322 !array 4323 ! ************************************************************************* 4324 4325 do i=1,n_list 4326 call sort_combination(combination_list(:,i),n_int) 4327 end do 4328 4329 end subroutine sort_combination_list