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