TABLE OF CONTENTS
ABINIT/compute_anharmonics [ Functions ]
NAME
compute_anharmonics
FUNCTION
Compute strain phonon coupling by finite differences Return the effective_potential with the third order
INPUTS
filenames(17) = path with all name files inp <type(multibinit_dtset_type)> = datatype with all the input variables comm=MPI communicator
OUTPUT
eff_pot<type(effective_potential_type)> = effective_potential datatype to be initialized
PARENTS
multibinit
CHILDREN
effective_potential_file_read,effective_potential_free
effective_potential_setelastic3rd,effective_potential_setelastic4th
effective_potential_setelasticdispcoupling
effective_potential_setstrainphononcoupling
effective_potential_writeabiinput,harmonics_terms_applysumrule,ifc_free
strain_free,strain_get,strain_init,wrtout,xmpi_bcast
SOURCE
32 #if defined HAVE_CONFIG_H 33 #include "config.h" 34 #endif 35 36 #include "abi_common.h" 37 38 subroutine compute_anharmonics(eff_pot,filenames,inp,comm) 39 40 use defs_basis 41 use m_errors 42 use m_profiling_abi 43 use m_xmpi 44 use m_io_tools, only : open_file 45 46 use m_ifc 47 use m_anharmonics_terms 48 use m_effective_potential 49 use m_effective_potential_file 50 use m_multibinit_dataset, only : multibinit_dtset_type 51 use m_strain 52 use m_fstrings, only : itoa,int2char4,ftoa 53 54 !This section has been created automatically by the script Abilint (TD). 55 !Do not modify the following lines by hand. 56 #undef ABI_FUNC 57 #define ABI_FUNC 'compute_anharmonics' 58 use interfaces_14_hidewrite 59 !End of the abilint section 60 61 implicit none 62 63 !Arguments ------------------------------------ 64 !scalars 65 integer, intent(in) :: comm 66 character(len=fnlen),intent(in) :: filenames(17) 67 type(effective_potential_type),target, intent(inout) :: eff_pot 68 type(multibinit_dtset_type),intent(in) :: inp 69 !arrays 70 71 !Local variables------------------------------- 72 !scalar 73 integer :: ia,ii,ierr,irpt,jj,kk,my_rank,natom 74 integer :: nfile,nrpt,nproc 75 real(dp) :: delta,delta1,delta2 76 character(len=500) :: message 77 character(len=fnlen):: name 78 logical :: files_availables = .True.,has_any_strain = .False. 79 logical :: has_all_strain = .True. 80 logical :: iam_master=.FALSE. 81 integer,parameter :: master=0 82 !arrays 83 integer :: have_strain(6) 84 real(dp) :: deformation(6,2),elastics3rd(6,6,6) 85 real(dp) :: elastics4th(6,6,6,6),rprimd_def(3,3) 86 type(strain_type) :: strain 87 type(ifc_type) :: phonon_strain(6) 88 logical, allocatable :: file_usable(:) 89 real(dp),allocatable :: elastic_displacement(:,:,:,:) 90 type(effective_potential_type),dimension(:),allocatable :: eff_pots 91 type(strain_type),dimension(:),allocatable :: effpot_strain 92 type(effective_potential_type),pointer :: ref_eff_pot 93 94 ! ************************************************************************* 95 96 write(message,'(a,(80a),a)') ch10,('=',ii=1,80),ch10 97 call wrtout(ab_out,message,'COLL') 98 call wrtout(std_out,message,'COLL') 99 100 write(message, '(a,a,a)' )' Compute the third order derivative by finite differences',ch10 101 call wrtout(std_out,message,'COLL') 102 call wrtout(ab_out,message,'COLL') 103 104 write(message, '(a,a,a)' )' The following files will be used :' 105 call wrtout(std_out,message,'COLL') 106 call wrtout(ab_out,message,'COLL') 107 108 !========================================== 109 !0)Initialisation of variables: 110 ! Set MPI local varibaless 111 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 112 iam_master = (my_rank == master) 113 114 !========================================== 115 !1) Get the list of files 116 nfile = 0 117 jj=6 118 do while (jj < 18) 119 if (filenames(jj)/="") then 120 if(jj==6) nfile = 0 121 write(message, '(a,a)' )' - ',trim(filenames(jj)) 122 call wrtout(std_out,message,'COLL') 123 call wrtout(ab_out,message,'COLL') 124 jj = jj + 1 125 nfile = nfile + 1 126 else 127 exit 128 end if 129 end do 130 131 if(nfile==0) then 132 write(message,'(a)') ' - No file found -' 133 call wrtout(ab_out,message,'COLL') 134 call wrtout(std_out,message,'COLL') 135 end if 136 137 write(message,'(a,(80a),a)') ch10,('-',ii=1,80),ch10 138 call wrtout(ab_out,message,'COLL') 139 call wrtout(std_out,message,'COLL') 140 141 !============================================ 142 !2) Read the effectives potential from files" 143 ! - store the reference effective potential 144 ! - Also get the strain 145 ! - perform some checks 146 ABI_DATATYPE_ALLOCATE(eff_pots,(nfile)) 147 ABI_DATATYPE_ALLOCATE(effpot_strain,(nfile)) 148 ABI_ALLOCATE(file_usable,(nfile)) 149 150 ref_eff_pot => eff_pot 151 file_usable(:) = .True. 152 153 ii = 1 ! Start at the index 1 154 jj = 6 ! Start at the index 6 155 do while (jj < 18) 156 if (filenames(jj)/="".and.filenames(jj)/="no") then 157 !Read and Intialisation of the effective potential type 158 call effective_potential_file_read(filenames(jj),eff_pots(ii),inp,comm) 159 !Eventualy print the xml file 160 ! if(inp%prt_model==-1.or.inp%prt_model>=3) then 161 ! call int2char4(ii,message) 162 ! name = 'structure_'//trim(itoa(ii-1))//'.xml' 163 ! call isfile(name,'new') 164 ! call effective_potential_writeXML(eff_pots(ii),1,filename=name) 165 ! end if 166 167 !Fill the eff_pots with the conresponding strain 168 call strain_get(effpot_strain(ii),rprim=eff_pot%crystal%rprimd,& 169 & rprim_def=eff_pots(ii)%crystal%rprimd) 170 171 jj = jj + 1; ii = ii + 1 172 173 write(message,'(a,(80a))') ch10,('-',ia=1,80) 174 call wrtout(ab_out,message,'COLL') 175 call wrtout(std_out,message,'COLL') 176 else 177 exit 178 end if 179 end do 180 181 !Do some checks 182 if(iam_master)then 183 do ii=1,size(eff_pots) 184 if (eff_pots(ii)%harmonics_terms%ifcs%nrpt/=ref_eff_pot%harmonics_terms%ifcs%nrpt) then 185 write(message,'(a,I0,a,a,a,a,a,I0,a,a,a,a)' )& 186 & 'the number of cell in reference (',ref_eff_pot%harmonics_terms%ifcs%nrpt,& 187 & ') is not equal to the ',ch10,'the number of cell in ',trim(filenames(ii+5)),& 188 & ' (',eff_pots(ii)%harmonics_terms%ifcs%nrpt,')',ch10,'this files cannot be used',ch10 189 MSG_WARNING(message) 190 file_usable(ii) = .False. 191 end if 192 if (eff_pots(ii)%crystal%natom/=ref_eff_pot%crystal%natom) then 193 write(message, '(a,I0,a,a,a,a,a,I0,a,a,a,a)' )& 194 & 'the number of atoms in reference (',ref_eff_pot%crystal%natom,') is not equal to the ',ch10,& 195 & 'the number of atoms in ',trim(filenames(ii+5)),' (',eff_pots(ii)%crystal%natom,')',ch10,& 196 & 'this files cannot be used',ch10 197 MSG_WARNING(message) 198 file_usable(ii) = .False. 199 end if 200 if (eff_pots(ii)%crystal%ntypat/=ref_eff_pot%crystal%ntypat) then 201 write(message, '(a,I0,a,a,a,a,a,I0,a,a,a,a)' )& 202 & 'the number of type of atoms in reference (',ref_eff_pot%crystal%ntypat,& 203 & ') is not equal to the ',& 204 & ch10,'the number of type of atoms in ',trim(filenames(ii+5)),& 205 & ' (',eff_pots(ii)%crystal%ntypat,')',& 206 & ch10,'this files can not be used',ch10 207 MSG_WARNING(message) 208 file_usable(ii) = .False. 209 end if 210 end do 211 end if 212 213 ! MPI BROADCAST 214 do ii=1,size(eff_pots) 215 call xmpi_bcast (file_usable(ii), master, comm, ierr) 216 end do 217 218 if (count((effpot_strain%name=="reference"))>1) then 219 write(message, '(2a)' )& 220 & ' There is several file corresponding to the reference ',ch10 221 MSG_BUG(message) 222 end if 223 224 have_strain = 0 225 226 write(message,'(a)') ' Strains available after reading the files:' 227 call wrtout(ab_out,message,'COLL') 228 call wrtout(std_out,message,'COLL') 229 do ii=1,size(eff_pots) 230 if(effpot_strain(ii)%name /= "".and.file_usable(ii)) then 231 write(message,'(a,a,a,I2,a,(ES10.2),a)')& 232 & ' A ',trim(effpot_strain(ii)%name),' strain in the direction ',& 233 & effpot_strain(ii)%direction,' with delta of ',effpot_strain(ii)%delta 234 has_any_strain = .True. 235 call wrtout(ab_out,message,'COLL') 236 call wrtout(std_out,message,'COLL') 237 end if 238 end do 239 240 241 if(nfile>1.and.has_any_strain) then 242 write(message,'(a,a,a)') ch10, ' ---analize in more details these files---',ch10 243 call wrtout(ab_out,message,'COLL') 244 call wrtout(std_out,message,'COLL') 245 else 246 write(message,'(a)') ' - No strain found -' 247 call wrtout(ab_out,message,'COLL') 248 call wrtout(std_out,message,'COLL') 249 write(message,'(a,(80a),a)') ch10,('-',ia=1,80),ch10 250 call wrtout(ab_out,message,'COLL') 251 call wrtout(std_out,message,'COLL') 252 end if 253 254 !First check the strain 255 do ii =1,6 256 jj = 0 257 jj = count(effpot_strain%direction==ii) 258 if(jj>2) then 259 write(message, '(a,I1,a)' )& 260 & ' There is several file corresponding to strain uniaxial in direction ',ii,ch10 261 MSG_ERROR(message) 262 else 263 name = 'uniaxial' 264 if(ii>=4) name = 'shear' 265 if (jj==1) then 266 write(message, '(a,a,a,I1,a,a)' )& 267 & ' WARNING: There is only one strain ',trim(name),' in direction ',ii,ch10,& 268 & ' the finate diferences will not be centering' 269 call wrtout(std_out,message,"COLL") 270 has_all_strain = .False. 271 have_strain(ii)=jj 272 else 273 if(jj==2)then 274 write(message, '(a,a,a,I1,a)' )& 275 & ' There is two files corresponding to strain ',trim(name),' in direction ',ii,ch10 276 call wrtout(ab_out,message,'COLL') 277 call wrtout(std_out,message,'COLL') 278 have_strain(ii)=jj 279 else 280 write(message, '(a,a,a,I1,a,a)' )& 281 & ' WARNING: There is no strain ',trim(name),' in direction ',ii,ch10 282 call wrtout(std_out,message,"COLL") 283 has_all_strain = .False. 284 if (inp%strcpling == 2) then 285 do kk = 1,2 286 delta = inp%delta_df 287 if (kk==1) delta = -1 * delta 288 call strain_init(strain,name=name,direction=ii,delta=delta) 289 rprimd_def = matmul(eff_pot%crystal%rprimd,strain%strain) 290 if(kk==1) then 291 write(message, '(a,a,a,a,a,I1,a,a,a,a)' )& 292 & ' if you want to get the correct structure, please run dfpt calculation with',ch10,& 293 & ' strain ',trim(name),' in the direction',ii,' with delta=',trim(ftoa(delta)),ch10,& 294 & ' The corresponding primitive vectors are:' 295 else 296 write(message, '(a,a,a,I1,a,a,a,a)' )& 297 & ' And a strain ',trim(name),' in the direction',ii,' with delta = ',& 298 & trim(ftoa(delta)),ch10,' The corresponding primitive vectors are:' 299 end if 300 call wrtout(ab_out,message,'COLL') 301 call wrtout(std_out,message,'COLL') 302 write(message,'(3(F20.10),a,3(F20.10),a,3(F20.10))')& 303 & rprimd_def(:,1),ch10, rprimd_def(:,2), ch10,rprimd_def(:,3) 304 call wrtout(ab_out,message,'COLL') 305 call wrtout(std_out,message,'COLL') 306 if(iam_master)then 307 call effective_potential_writeAbiInput(eff_pot,strain=strain) 308 end if 309 call strain_free(strain) 310 end do 311 end if 312 end if 313 end if 314 end if 315 end do 316 317 ! check if strain exist 318 if(all(have_strain==0).and.inp%strcpling /= 2) then 319 write(message, '(6a)' )& 320 & ' WARNING: There is no file corresponding to strain',& 321 & ' to compute 3rd order derivatives.',ch10,& 322 & ' In this case the 3rd order derivatives are not set',ch10,& 323 & ' (add files or set strcpling to 0)' 324 call wrtout(std_out,message,"COLL") 325 end if 326 327 ! check if the existing strains have opposite deformation 328 deformation = zero 329 do ii=1,6 330 if(have_strain(ii)/=0) then 331 ia = 1 332 do jj=1,size(eff_pots) 333 if (effpot_strain(jj)%direction==ii)then 334 deformation(ii,ia) = effpot_strain(jj)%delta 335 ia = ia + 1 336 end if 337 end do 338 if (have_strain(ii)==2) then 339 delta1 = deformation(ii,1) 340 delta2 = deformation(ii,2) 341 if (delta1+delta2 > tol15) then 342 write(message, '(a,I1,a,a)' )& 343 & ' The deformations for strain ',ii,& 344 & ' are not the opposite',ch10 345 MSG_ERROR(message) 346 end if 347 end if 348 end if 349 end do 350 351 write(message,'(a,(80a))') ch10,('-',ia=1,80) 352 call wrtout(ab_out,message,'COLL') 353 call wrtout(std_out,message,'COLL') 354 355 write(message,'(a,a)') ch10, ' After analyzing, the strains available are:' 356 call wrtout(ab_out,message,'COLL') 357 call wrtout(std_out,message,'COLL') 358 if(has_any_strain) then 359 do ii=1,6 360 if(have_strain(ii)/=0) then 361 do jj=1,size(eff_pots) 362 if (effpot_strain(jj)%direction==ii)then 363 write(message,'(a,a,a,I2,a,(ES10.2),a)')& 364 & ' A ',trim(effpot_strain(jj)%name),' strain in the direction ',& 365 & effpot_strain(jj)%direction,' with delta of ',effpot_strain(jj)%delta 366 call wrtout(ab_out,message,'COLL') 367 call wrtout(std_out,message,'COLL') 368 end if 369 end do 370 else 371 files_availables = .False. 372 end if 373 end do 374 else 375 files_availables = .False. 376 write(message,'(a)') ' - No strain available -' 377 call wrtout(ab_out,message,'COLL') 378 call wrtout(std_out,message,'COLL') 379 end if 380 write(message,'(a,(80a))') ch10,('-',ia=1,80) 381 call wrtout(ab_out,message,'COLL') 382 call wrtout(std_out,message,'COLL') 383 384 if(has_all_strain) then 385 write(message,'(3a)') ch10, ' The computation of the third order derivative ',& 386 & 'is possible' 387 else 388 if (inp%strcpling /= 2) then 389 if(ref_eff_pot%has_anharmonicsTerms)then 390 write(message,'(10a)') ch10, ' The computation of the third order derivative ',& 391 & 'is not possible',ch10,' somes files are missing please use strcpling 2 to generate',& 392 & ' inputs files',ch10,' usable by abinit. The third order derivatives present in ',& 393 & trim(filenames(3)),' will be used' 394 else 395 write(message,'(9a)') ch10, ' The computation of the third order derivative ',& 396 & 'is not possible',ch10,' somes files are missing please use strcpling 2 to generate',& 397 & ' inputs files',ch10,' usable by abinit. The third order derivative will not be set in',& 398 & ' the XML file' 399 end if 400 else 401 if(ref_eff_pot%has_anharmonicsTerms)then 402 write(message,'(10a)') ch10, ' The computation of the third order derivative ',& 403 & 'is not possible',ch10,' somes files are missing, the input files usable by abinit have been',& 404 & ' generate.',ch10,' The third order derivatives present in ',trim(filenames(3)),' will be used' 405 else 406 write(message,'(8a)') ch10, ' The computation of the third order derivative ',& 407 & 'is not possible',ch10,' somes files are missing, the input files usable by abinit have been',& 408 & ' generate.',ch10,' The third order derivatives will be not set in the XML file' 409 end if 410 end if 411 call wrtout(ab_out,message,'COLL') 412 call wrtout(std_out,message,'COLL') 413 end if 414 415 !================================================ 416 !3) Compute finate differences 417 if(has_all_strain) then 418 419 ! Allocation of array and set some values 420 nrpt = ref_eff_pot%harmonics_terms%ifcs%nrpt 421 natom = ref_eff_pot%crystal%natom 422 ABI_ALLOCATE(elastic_displacement,(6,6,3,natom)) 423 424 elastics3rd = zero 425 elastics4th = zero 426 427 do ii=1,6 428 if(have_strain(ii)/=0) then 429 ! We want the find the index of the perturbation ii in eff_pots(ii) 430 ! And store in delta1 and delta2 431 delta1 = zero 432 delta2 = zero 433 do jj=1,size(eff_pots) 434 if (effpot_strain(jj)%direction==ii.and.(effpot_strain(jj)%direction/=0))then 435 if (abs(delta1)<tol16) then 436 delta1 = jj 437 else 438 delta2 = jj 439 end if 440 end if 441 end do 442 if (abs(delta1)>tol16.and.abs(delta1)>tol16)then 443 ! check if delta1 < delta2, in this case, inverse delta1 and delta2 444 if (effpot_strain(int(delta1))%delta < effpot_strain(int(delta2))%delta) then 445 delta = delta1 446 delta1 = delta2 447 delta2 = delta 448 end if 449 ! Compute strain phonon-coupling 450 phonon_strain(ii)%nrpt = nrpt 451 ABI_ALLOCATE(phonon_strain(ii)%atmfrc,(3,natom,3,natom,nrpt)) 452 ABI_ALLOCATE(phonon_strain(ii)%cell,(3,nrpt)) 453 phonon_strain(ii)%atmfrc = zero 454 phonon_strain(ii)%cell = eff_pots(int(delta1))%harmonics_terms%ifcs%cell 455 456 do irpt=1,phonon_strain(ii)%nrpt 457 phonon_strain(ii)%atmfrc(:,:,:,:,irpt) =& 458 & (eff_pots(int(delta1))%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt)& 459 & - eff_pots(int(delta2))%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt)) / & 460 & (2 * abs(effpot_strain(int(delta1))%delta)) 461 end do 462 463 if(inp%asr >= 0) then 464 ! Impose sum rule 465 call harmonics_terms_applySumRule(inp%asr,phonon_strain(ii),& 466 & eff_pot%crystal%natom) 467 end if 468 469 ! Compute elastic constants 470 elastics3rd(ii,:,:) = (eff_pots(int(delta1))%harmonics_terms%elastic_constants(:,:)& 471 & - eff_pots(int(delta2))%harmonics_terms%elastic_constants(:,:)) / & 472 & (2 * abs(effpot_strain(int(delta1))%delta)) 473 474 ! Compute elastic-displacement coupling 475 elastic_displacement(ii,:,:,:)=(eff_pots(int(delta1))%harmonics_terms%strain_coupling(:,:,:)& 476 & - eff_pots(int(delta2))%harmonics_terms%strain_coupling(:,:,:)) / & 477 & (2 * abs(effpot_strain(int(delta1))%delta)) 478 479 ! Compute elastic constants 480 elastics4th(ii,ii,:,:) = (eff_pots(int(delta1))%harmonics_terms%elastic_constants(:,:)& 481 & - 2*ref_eff_pot%harmonics_terms%elastic_constants(:,:)& 482 & + eff_pots(int(delta2))%harmonics_terms%elastic_constants(:,:)) / & 483 & (abs(effpot_strain(int(delta1))%delta)**2) 484 end if 485 486 end if 487 end do 488 489 ! Set all the values in the effective potential type 490 call effective_potential_setStrainPhononCoupling(eff_pot,natom,phonon_strain) 491 call effective_potential_setElastic3rd(eff_pot,elastics3rd) 492 call effective_potential_setElastic4th(eff_pot,elastics4th) 493 call effective_potential_setElasticDispCoupling(eff_pot,natom,elastic_displacement) 494 495 496 ! Free the phonon-strain coupling array 497 do ii = 1,6 498 call ifc_free(phonon_strain(ii)) 499 end do 500 ABI_DEALLOCATE(elastic_displacement) 501 502 write(message,'(4a)') ch10, ' The computation of the 3rd order elastics constants, ',ch10,& 503 & ' the phonon-strain coupling and the elastic-displacement coupling is done' 504 call wrtout(ab_out,message,'COLL') 505 call wrtout(std_out,message,'COLL') 506 507 end if 508 509 510 !=============================================== 511 !4) Free the array of effective potential 512 513 do jj=1,nfile 514 ! Free the effective potential type 515 call effective_potential_free(eff_pots(jj)) 516 end do 517 518 ABI_DATATYPE_DEALLOCATE(effpot_strain) 519 ABI_DATATYPE_DEALLOCATE(eff_pots) 520 ABI_DEALLOCATE(file_usable) 521 522 523 write(message,'(a,a,a,(80a))') ch10,('=',ii=1,80),ch10 524 call wrtout(ab_out,message,'COLL') 525 call wrtout(std_out,message,'COLL') 526 527 end subroutine compute_anharmonics