TABLE OF CONTENTS
- ABINIT/m_effective_potential
- m_effective_potential/effective_potential_checkDEV
- m_effective_potential/effective_potential_copy
- m_effective_potential/effective_potential_init
- m_effective_potential/effective_potential_initmpi
- m_effective_potential/effective_potential_setConfinement
- m_effective_potential/effective_potential_setSupercell
- m_effective_potential/effective_potential_type
- m_effective_potential/effective_potential_writeAbiInput
- m_effective_potential/effective_potential_writeNETCDF
- m_effective_potential/effective_potential_writeXML
- m_effective_potential/equal
ABINIT/m_effective_potential [ Modules ]
NAME
m_effective_potential
FUNCTION
Module for the effective potential Container type is defined, and destruction, print subroutines Contain also routine to evaluate the energy,forces and stresses
COPYRIGHT
Copyright (C) 2010-2024 ABINIT group (AM) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 module m_effective_potential 27 28 use defs_basis 29 use defs_abitypes 30 use m_errors 31 use m_abicore 32 use m_strain 33 use m_ifc 34 use m_supercell 35 use m_phonons 36 use m_ddb 37 use m_polynomial_conf 38 use m_polynomial_coeff 39 use m_anharmonics_terms 40 use m_harmonics_terms 41 use m_xmpi 42 use m_ewald 43 use m_nctk 44 #if defined HAVE_NETCDF 45 use netcdf 46 #endif 47 #if defined DEV_MS_SCALEUP 48 use scup_global, only : global_calculate_energy, global_calculate_forces 49 #endif 50 51 use m_fstrings, only : replace, ftoa, itoa 52 use m_io_tools, only : open_file, get_unit 53 use m_dtfil, only : isfile 54 use m_symtk, only : matr3inv 55 use m_effpot_mpi, only : effpot_mpi_init,effpot_mpi_type,effpot_mpi_free 56 use m_abihist, only : abihist 57 use m_geometry, only : gred2fcart,fcart2gred, xcart2xred, xred2xcart, metric 58 use m_crystal, only : crystal_t, crystal_init 59 !use m_anaddb_dataset, only : anaddb_dataset_type, anaddb_dtset_free, outvars_anaddb, invars9 60 61 implicit none 62 63 public :: effective_potential_distributeResidualForces 64 public :: effective_potential_evaluate 65 public :: effective_potential_free 66 public :: effective_potential_freeCoeffs 67 public :: effective_potential_freempi 68 public :: effective_potential_generateDipDip 69 public :: effective_potential_getDisp 70 public :: effective_potential_init 71 public :: effective_potential_initmpi 72 public :: effective_potential_copy 73 public :: effective_potential_print 74 public :: effective_potential_printSupercell 75 public :: effective_potential_setCoeffs 76 public :: effective_potential_setConfinement 77 public :: effective_potential_setElastic3rd 78 public :: effective_potential_setElastic4th 79 public :: effective_potential_setElasticDispCoupling 80 public :: effective_potential_setStrainPhononCoupling 81 public :: effective_potential_setSupercell 82 public :: effective_potential_writeAbiInput 83 public :: effective_potential_writeXML 84 public :: effective_potential_writeAnhHead 85 !AM_EXPERIMENTAL 86 public :: effective_potential_computeGradient 87 ! public :: effective_potential_effpot2ddb 88 ! public :: effective_potential_printPDOS 89 public :: effective_potential_checkDEV 90 public :: effective_potential_writeNETCDF 91 public :: OPERATOR(==) 92 !AM_EXPERIMENTAL
m_effective_potential/effective_potential_checkDEV [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_checkDEV
FUNCTION
Routine for develloper Check by finite differences the equations in effective_potential_evaluate need to provide HIST file, so you need to activate the fit_process or bound_process to activate the reading of the HIST
INPUTS
eff_pot<type(effective_potential)> = effective potential hist<type(abihist)> = The history of the MD natom = number of atom ntime = number of time in the hist
OUTPUT
SOURCE
3411 subroutine effective_potential_checkDEV(eff_pot,hist,natom,ntime) 3412 3413 !Arguments ------------------------------------ 3414 !scalars 3415 integer, intent(in) :: natom,ntime 3416 !arrays 3417 type(effective_potential_type),intent(in) :: eff_pot 3418 type(abihist),intent(in) :: hist 3419 !Local variables------------------------------- 3420 !scalar 3421 integer :: ii,jj,ia,mu,npt,istep 3422 ! integer :: ifirst 3423 real(dp):: energy,delt,delta,ucvol 3424 !arrays 3425 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),mat_def(3,3),identity(3,3) 3426 real(dp):: fcart(3,natom),gred(3,natom),strten(6),rprimd(3,3) 3427 real(dp):: rprimd_def(3,3),rprimd_ref(3,3),deltalist(5) 3428 real(dp):: disp(3,natom),disp_red(3,natom),strain(6),du_delta(6,3,natom),diff(5) 3429 real(dp),allocatable :: xred(:,:) 3430 integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/) 3431 character(len=500) :: msg 3432 3433 ! ************************************************************************* 3434 3435 !Do some checks 3436 if(ntime /= hist%mxhist)then 3437 write(msg,'(a)')'ntime is not correct' 3438 ABI_BUG(msg) 3439 end if 3440 3441 if(natom /= size(hist%xred,2)) then 3442 write(msg,'(a)')'natom is not correct' 3443 ABI_BUG(msg) 3444 end if 3445 3446 3447 ABI_MALLOC(xred,(3,natom)) 3448 xred = zero 3449 3450 !option 1 => set the reference for the test 3451 ! call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,& 3452 !& eff_pot%supercell%xcart,xred) 3453 ! rprimd = eff_pot%supercell%rprimd 3454 3455 !option 2 => set a specific step for the test 3456 istep = 4 3457 xred = hist%xred(:,:,istep) 3458 rprimd = hist%rprimd(:,:,istep) 3459 3460 rprimd_ref = eff_pot%supercell%rprimd 3461 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 3462 3463 npt=5 3464 delta = 0.001 3465 deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/) 3466 strain = zero 3467 3468 do ia=1,natom 3469 do mu=1,3 3470 write(std_out,*) "atm: ",ia," dir: ",mu 3471 do ii=1,npt 3472 delt = deltalist(ii) 3473 3474 ! Get the initial displacement 3475 call effective_potential_getDisp(disp,du_delta,natom,rprimd,& 3476 & eff_pot%supercell%rprimd,1,xred_hist=xred,& 3477 & xcart_ref=eff_pot%supercell%xcart,& 3478 & compute_displacement = .true.,compute_duDelta = .true.) 3479 3480 ! Add the delta 3481 call xcart2xred(natom, rprimd, disp, disp_red) 3482 disp_red(mu,ia) = disp_red(mu,ia) + delt 3483 call xred2xcart(natom, rprimd, disp, disp_red) 3484 3485 call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,& 3486 & xred=xred,du_delta=du_delta,& 3487 & displacement=disp,compute_anharmonic=.true.,verbose=.false.) 3488 diff(ii) = energy 3489 3490 end do 3491 3492 ! Get the initial displacement 3493 call effective_potential_getDisp(disp,du_delta,natom,rprimd,& 3494 & eff_pot%supercell%rprimd,1,xred_hist=xred,& 3495 & xcart_ref=eff_pot%supercell%xcart,& 3496 & compute_displacement = .true.,compute_duDelta = .true.) 3497 3498 call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,& 3499 & xred=xred,du_delta=du_delta,& 3500 & displacement=disp,compute_anharmonic=.true.,verbose=.false.) 3501 3502 write(std_out,*) "Analyti:",gred(mu,ia) 3503 write(std_out,*) "FD :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) 3504 write(std_out,*) "Diff(%):",abs(100*(gred(mu,ia)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))& 3505 & / (12*delta) )) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) )) 3506 3507 end do 3508 end do 3509 3510 3511 ! Fill the identity matrix 3512 identity = zero 3513 forall(ii=1:3)identity(ii,ii)=1 3514 3515 npt=5 3516 delta = 0.0005 3517 deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/) 3518 3519 do jj=1,6 3520 write(std_out,*) "strain ",jj 3521 do ii=1,npt 3522 strain = zero 3523 delt = deltalist(ii) 3524 mat_def = zero 3525 strain(jj) = strain(jj) + delt 3526 3527 mat_def(alpha(jj),beta(jj)) = mat_def(alpha(jj),beta(jj)) + half * strain(jj) 3528 mat_def(beta(jj),alpha(jj)) = mat_def(beta(jj),alpha(jj)) + half * strain(jj) 3529 3530 rprimd_def = matmul(identity(:,:)+mat_def(:,:),rprimd) 3531 3532 ! The two options should give the same result 3533 ! Option 1 => compute the disps and provide them to evaluate 3534 ! call effective_potential_getDisp(disp,du_delta,natom,rprimd_def,& 3535 ! & rprimd_ref,1,xred_hist=xred,& 3536 ! & xcart_ref=eff_pot%supercell%xcart,& 3537 ! & compute_displacement = .true.,compute_duDelta = .true.) 3538 3539 ! call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd_def,& 3540 ! & xred=xred,du_delta=du_delta,& 3541 ! & displacement=disp,strain=strain,& 3542 ! & compute_anharmonic=.true.,verbose=.false.) 3543 3544 ! Option 2 => compute the disps within evaluate 3545 call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd_def,& 3546 & xred=xred,compute_anharmonic=.true.,verbose=.false.) 3547 3548 3549 3550 diff(ii) = energy 3551 3552 end do 3553 3554 ! The two options should give the same result 3555 ! Option 1 => compute the disps and provide them to evaluate 3556 ! call effective_potential_getDisp(disp,du_delta,natom,rprimd,& 3557 ! & rprimd_ref,1,xred_hist=xred,& 3558 ! & xcart_ref=eff_pot%supercell%xcart,& 3559 ! & compute_displacement = .true.,compute_duDelta = .true.) 3560 3561 ! call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,& 3562 ! & xred=xred,du_delta=du_delta,& 3563 ! & displacement=disp,& 3564 ! & compute_anharmonic=.true.,verbose=.false.) 3565 3566 ! Option 2 => compute the disps within evaluate 3567 call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,& 3568 & xred=xred,compute_anharmonic=.true.,verbose=.false.) 3569 3570 write(std_out,*) "Analyti:",strten(jj) 3571 write(std_out,*) "FD :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol 3572 write(std_out,*) "Diff(%):",abs(100*(strten(jj)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))& 3573 & / (12*delta) / ucvol)) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol)) 3574 3575 end do 3576 3577 end subroutine effective_potential_checkDEV
m_effective_potential/effective_potential_copy [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_copy
FUNCTION
Copy one effective potential to another i.e. initialize the output effective potentiale eff_pot_out with the values the input effective potential eff_pot_in
INPUTS
eff_pot_in<type(effective_potential_type)> = effective_potential datatype to be initialized
OUTPUT
eff_pot_out<type(effective_potential_type)> = effective_potential datatype to be initialized
SOURCE
453 subroutine effective_potential_copy(eff_pot_out,eff_pot_in,comm) 454 455 !Arguments ------------------------------------ 456 !scalars 457 type(effective_potential_type),intent(out) :: eff_pot_out 458 type(effective_potential_type),intent(in) :: eff_pot_in 459 integer,intent(in) :: comm 460 ! *********************************************************************** 461 462 call effective_potential_init(eff_pot_in%crystal,eff_pot_out,eff_pot_in%energy,eff_pot_in%harmonics_terms%ifcs,& 463 & eff_pot_in%anharmonics_terms%ncoeff,eff_pot_in%harmonics_terms%nqpt,comm,& 464 & coeffs=eff_pot_in%anharmonics_terms%coefficients,dynmat=eff_pot_in%harmonics_terms%dynmat,& 465 & elastic_constants=eff_pot_in%harmonics_terms%elastic_constants,& 466 & elastic3rd=eff_pot_in%anharmonics_terms%elastic3rd,& 467 & elastic_displacement=eff_pot_in%anharmonics_terms%elastic_displacement,& 468 & epsilon_inf=eff_pot_in%harmonics_terms%epsilon_inf,& 469 & fcart=eff_pot_in%fcart,strain_coupling=eff_pot_in%harmonics_terms%strain_coupling,& 470 & strten=eff_pot_in%strten,name=eff_pot_in%name,& 471 & phonon_strain=eff_pot_in%anharmonics_terms%phonon_strain,phfrq=eff_pot_in%harmonics_terms%phfrq,& 472 & qpoints=eff_pot_in%harmonics_terms%qpoints,has_anharmonicsTerms=eff_pot_in%has_anharmonicsTerms,& 473 & supercell=eff_pot_in%supercell,& 474 & zeff=eff_pot_in%harmonics_terms%zeff) 475 476 477 end subroutine effective_potential_copy
m_effective_potential/effective_potential_init [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_init
FUNCTION
Initialize effective_potential datatype
INPUTS
crytal<type(crystal_t)> = datatype with all the information for the crystal energy = energy of the reference structure ifcs <type(ifc_type)> = ifc type with cell,ewald short and total range of the ifcs ncoeff = number of coefficients in the polynomial nqpt = number of qpoints comm = mpi comunicator coeffs(ncoeff)<type(polynomial_coeff_type)> = optional,list of coefficients for the anharmonic part dynmat(2,3,natom,3,natom,nqpt) = optional,dynamical matrix for each qpoints epsilon_inf(3,3) = optional,dielectric tensor elastic_constants(6,6) = optional,elastic constants tensor elastic3rd(6,6,6) = optional,3 order derivatives with respect to to 3 strain elastic_displacement(6,6,3,natom)=optional, 3 order derivatives with respect to 2 strain and 1 Atom disp fcart(3,natom) = optional,optional,initial fcart in the structure strain_coupling(6,natom,3) = optional, internal strain coupling parameters strten(6) = optional,optional,initial strain in the structure name = optional, name of the structure phonon_strain(6) = optional,ifc type for the phonon-strain coupling (should be in anharmonics_terms) phfrq(3*natom,nqpt) = optional,phonons frequencies for each q points in Hartree/cm qpoints(3,nqpt) = optional,list of qpoints wavevectors has_anharmonicsTerms = optional, (default false) flag to set the anharmonics terms supercell<type(supercell_type)> = optional, supercell type to define zeff(3,natom) = optional,effective charges
OUTPUT
eff_pot<type(effective_potential_type)> = effective_potential datatype to be initialized
SOURCE
195 subroutine effective_potential_init(crystal,eff_pot,energy,ifcs,ncoeff,nqpt,comm,& 196 & coeffs,dynmat,elastic_constants,elastic3rd,& 197 & elastic_displacement,epsilon_inf,& 198 & fcart,strain_coupling,strten,name,phonon_strain,& 199 & polynomial_conf,phfrq,qpoints,has_anharmonicsTerms,& 200 & supercell,zeff) 201 202 !Arguments ------------------------------------ 203 !scalars 204 integer,intent(in) :: comm,ncoeff,nqpt 205 real(dp),intent(in):: energy 206 logical,optional,intent(in) :: has_anharmonicsTerms 207 !arrays 208 type(crystal_t),intent(in) :: crystal 209 type(effective_potential_type), intent(out) :: eff_pot 210 type(ifc_type),intent(in) :: ifcs 211 character(len=fnlen), optional,intent(in) :: name 212 real(dp),optional,intent(in) :: epsilon_inf(3,3),elastic_constants(6,6) 213 real(dp),optional,intent(in) :: dynmat(:,:,:,:,:,:),qpoints(:,:),phfrq(:,:) 214 real(dp),optional,intent(in) :: strain_coupling(:,:,:),zeff(:,:,:) 215 type(supercell_type),optional,intent(in) :: supercell 216 type(ifc_type),optional,intent(in) :: phonon_strain(6) 217 real(dp),optional,intent(in) :: elastic3rd(6,6,6) 218 real(dp),optional,intent(in) :: elastic_displacement(6,6,3,crystal%natom) 219 type(polynomial_coeff_type),optional,intent(in) :: coeffs(:) 220 real(dp),optional,intent(in) :: strten(6) 221 real(dp),optional,intent(in) :: fcart(3,crystal%natom) 222 type(polynomial_conf_type),optional,intent(in) :: polynomial_conf 223 224 !Local variables------------------------------- 225 !scalar 226 character(len=500) :: msg 227 !arrays 228 229 ! ************************************************************************* 230 231 !1-Free the effective potential before filling it 232 call effective_potential_free(eff_pot) 233 234 if (present(name)) then 235 eff_pot%name = name 236 else 237 eff_pot%name = '' 238 end if 239 240 !2-Perform some checks 241 if (crystal%natom < 1) then 242 write(msg, '(a,a,a,i10,a)' )& 243 & 'The cell must have at least one atom.',ch10,& 244 & 'The number of atom is ',crystal%natom,'.' 245 ABI_BUG(msg) 246 end if 247 248 if (crystal%ntypat < 1) then 249 write(msg, '(a,a,a,i10,a)' )& 250 & 'The cell must have at least one type of atom.',ch10,& 251 & 'The number of type of atom is ',crystal%ntypat,'.' 252 ABI_BUG(msg) 253 end if 254 255 !3-Fill energy of the crystal (hartree) 256 eff_pot%energy = energy 257 258 !1-Fill the crystal 259 !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here 260 call crystal_init(crystal%amu,eff_pot%crystal,crystal%space_group,crystal%natom,& 261 & crystal%npsp,crystal%ntypat,crystal%nsym,crystal%rprimd,& 262 & crystal%typat,crystal%xred,crystal%zion,crystal%znucl,& 263 & crystal%timrev,.FALSE.,.FALSE.,crystal%title,& 264 & symrel=crystal%symrel,tnons=crystal%tnons,symafm=crystal%symafm) 265 266 !4-Fill harmonic part 267 call harmonics_terms_init(eff_pot%harmonics_terms,ifcs,crystal%natom,ifcs%nrpt) 268 269 !5-Init the anharmonics_terms to set the flag to false 270 call anharmonics_terms_init(eff_pot%anharmonics_terms,crystal%natom,ncoeff) 271 272 !5-Fill optional inputs 273 ABI_MALLOC(eff_pot%fcart,(3,eff_pot%crystal%natom)) 274 eff_pot%fcart = zero 275 if(present(fcart))then 276 eff_pot%fcart = fcart 277 end if 278 279 if(present(elastic_constants))then 280 eff_pot%harmonics_terms%elastic_constants(:,:) = elastic_constants 281 end if 282 283 if(present(epsilon_inf))then 284 eff_pot%harmonics_terms%epsilon_inf(:,:) = epsilon_inf(:,:) 285 end if 286 287 if(present(dynmat).and.present(qpoints).and.present(phfrq))then 288 call harmonics_terms_setDynmat(dynmat,eff_pot%harmonics_terms,crystal%natom,nqpt,phfrq,qpoints) 289 end if 290 291 if(present(strain_coupling))then 292 call harmonics_terms_setInternalStrain(eff_pot%harmonics_terms,crystal%natom,strain_coupling) 293 end if 294 295 if(present(zeff))then 296 call harmonics_terms_setEffectiveCharges(eff_pot%harmonics_terms,crystal%natom,zeff) 297 end if 298 299 eff_pot%strten = zero 300 if(present(strten))then 301 eff_pot%strten(:) = strten(:) 302 end if 303 304 !Set the flag for the strain coupling 305 if(present(has_anharmonicsTerms)) then 306 eff_pot%has_anharmonicsTerms = has_anharmonicsTerms 307 else 308 eff_pot%has_anharmonicsTerms = .false. 309 end if 310 311 !Allocation of phonon strain coupling array (3rd order) 312 if(present(phonon_strain).and.has_anharmonicsTerms) then 313 call anharmonics_terms_setStrainPhononCoupling(eff_pot%anharmonics_terms,crystal%natom,& 314 & phonon_strain) 315 end if 316 317 !Set the 3rd order elastic tensor 318 if(present(elastic3rd).and.has_anharmonicsTerms)then 319 call anharmonics_terms_setElastic3rd(eff_pot%anharmonics_terms,elastic3rd) 320 end if 321 322 !TODO Comment Marcus: Uncomment and implement below if you want to use the elastic4th implementation, which is there, but 323 ! not used anywere 324 ! Below I just ensured that the has_elastic4th variable is set to FALSE as it was causing problems on some builders 325 ! in the test farm 326 ! Set the 4th order elastic tensor 327 ! if(present(elastic4th).and.has_anharmonicsTerms)then 328 ! call anharmonics_terms_setElastic3rd(eff_pot%anharmonics_terms,elastic4th) 329 ! end if 330 331 !MS Ensure has_elastic4th is .FALSE. 332 ! eff_pot%anharmonics_terms%has_elastic4th=.FALSE. 333 334 335 !Allocation of 3rd order with respecto to 2 strain and 1 atomic displacement 336 if(present(elastic_displacement).and.has_anharmonicsTerms)then 337 call anharmonics_terms_setElasticDispCoupling(eff_pot%anharmonics_terms,crystal%natom,& 338 & elastic_displacement) 339 end if 340 341 !Allocation of the coefficients 342 if(present(coeffs))then 343 if(ncoeff /= size(coeffs))then 344 ABI_BUG('ncoeff has not the same size than coeffs array') 345 end if 346 call effective_potential_setCoeffs(coeffs,eff_pot,ncoeff) 347 end if 348 349 if(present(supercell))then 350 call effective_potential_setSupercell(eff_pot,comm,supercell=supercell) 351 else 352 ! call init_supercell(eff_pot%crystal%natom, (/1,0,0, 0,1,0, 0,0,1/), eff_pot%crystal%rprimd,& 353 !& eff_pot%crystal%typat, eff_pot%crystal%xcart, eff_pot%crystal%znucl, eff_pot%supercell) 354 call effective_potential_setSupercell(eff_pot,comm,ncell=(/1,1,1/)) 355 end if 356 357 !Set the confinement potential 358 if(present(polynomial_conf)) then 359 call effective_potential_setConfinement(polynomial_conf%cutoff_disp,polynomial_conf%cutoff_strain,& 360 & eff_pot,polynomial_conf%factor_disp,& 361 & polynomial_conf%factor_strain,polynomial_conf%ndisp,& 362 & polynomial_conf%power_disp,polynomial_conf%power_strain,& 363 & polynomial_conf%need_confinement) 364 end if 365 366 end subroutine effective_potential_init
m_effective_potential/effective_potential_initmpi [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_initmpi
FUNCTION
Initializes the mpi information for parallelism over supercell. Only the parallelisation over cell is done here. The parallelisation over cell and coeff is disable for now (experimental)
INPUTS
eff_pot<type(effective_potential_type)> = datatype for the effective potential comm = MPI communicator
OUTPUT
This is for the parallelisation over the supercell eff_pot%mpi_ifc%me_supercell = Index of my processor in the comm. over one cell eff_pot%mpi_ifc%my_ncell = Number of cell treated by current proc eff_pot%mpi_ifc%my_cells(:) = Number of the cells in the supercell treat by this CPU eff_pot%mpi_ifc%my_index_cells(:,:) = indexes of the cells in the supercell treat by this CPU
SOURCE
392 subroutine effective_potential_initmpi(eff_pot,comm) 393 394 !Arguments ------------------------------------ 395 !scalars 396 type(effective_potential_type),intent(inout) :: eff_pot 397 integer,intent(in) :: comm 398 !arrays 399 400 !Local variables------------------------------- 401 !scalars 402 integer :: ndiv 403 integer :: ncell 404 !array 405 integer :: cell_number(3) 406 !character(len=500) :: msg 407 ! *********************************************************************** 408 409 !Set the number of cell in the supercell 410 cell_number(1) = eff_pot%supercell%rlatt(1,1) 411 cell_number(2) = eff_pot%supercell%rlatt(2,2) 412 cell_number(3) = eff_pot%supercell%rlatt(3,3) 413 ncell = product(cell_number(:)) 414 415 !Do some checks 416 if (any(cell_number <= 0).or.ncell<=0) then 417 ABI_ERROR('No supercell found for setting') 418 end if 419 420 !First mpi_ifc 421 ndiv = 1 422 call effpot_mpi_free(eff_pot%mpi_ifc) 423 call effpot_mpi_init(eff_pot%harmonics_terms%ifcs%cell,cell_number,eff_pot%mpi_ifc,& 424 & eff_pot%crystal%natom,ndiv,eff_pot%harmonics_terms%ifcs%nrpt,comm) 425 426 !Second mpi_coeff 427 ndiv = 1 428 call effpot_mpi_free(eff_pot%mpi_coeff) 429 call effpot_mpi_init(eff_pot%harmonics_terms%ifcs%cell,cell_number,eff_pot%mpi_coeff,& 430 & eff_pot%crystal%natom,ndiv,eff_pot%harmonics_terms%ifcs%nrpt,comm) 431 432 end subroutine effective_potential_initmpi
m_effective_potential/effective_potential_setConfinement [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_setConfinement
FUNCTION
Set the confinement in the effective_potential datatype
INPUTS
cutoff_disp(6) = Cutoff array for the strain cutoff_strain(ndisp) = Cutoff array for the atomic displacement factor_disp = Factor to appy to the polynomial term of the confinement (displacement) factor_strain = Factor to appy to the polynomial term of the confinement (strain) ndisp = Number of displacement (atoms) for the cut off power_disp = Power of the polynome related to the displacement power_strain = Power of the polynome related to the strain need_confinement = optional,Logical related to the necessity of the confinement
OUTPUT
eff_pot<type(effective_potential_type)> = datatype for effective potential
SOURCE
1307 subroutine effective_potential_setConfinement(cutoff_disp,cutoff_strain,eff_pot,factor_disp,& 1308 & factor_strain,ndisp,power_disp,power_strain,& 1309 & need_confinement) 1310 1311 !Arguments ------------------------------------ 1312 !scalars 1313 integer, intent(in) :: power_disp,power_strain,ndisp 1314 real(dp),intent(in) :: factor_disp,factor_strain 1315 logical,optional,intent(in) :: need_confinement 1316 !arrays 1317 real(dp),intent(in) :: cutoff_disp(ndisp),cutoff_strain(6) 1318 type(effective_potential_type),intent(inout) :: eff_pot 1319 !Local variables------------------------------- 1320 !scalar 1321 logical :: need_confinement_tmp 1322 !arrays 1323 !character(len=500) :: msg 1324 1325 ! ************************************************************************* 1326 1327 !Checks 1328 if (ndisp <= 0) then 1329 ABI_ERROR('ndisp can not be inferior or equal to zero') 1330 end if 1331 1332 !First free the type 1333 call polynomial_conf_free(eff_pot%confinement) 1334 1335 need_confinement_tmp = .FALSE. 1336 if (present(need_confinement)) need_confinement_tmp = need_confinement 1337 1338 call polynomial_conf_init(cutoff_disp,cutoff_strain,factor_disp,factor_strain,ndisp,& 1339 & eff_pot%confinement,power_disp,power_strain,& 1340 & need_confinement=need_confinement_tmp) 1341 1342 1343 end subroutine effective_potential_setConfinement
m_effective_potential/effective_potential_setSupercell [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_setSupercell
FUNCTION
Set the supercell type in the effective_potential type
INPUTS
eff_pot<type(effective_potential_type)> = effective_potential datatype comm = MPI communicator ncell(3) = optional, size of the supercell supercell = optional, supercell type to set to eff_pot
OUTPUT
eff_pot<type(effective_potential_type)> = effective_potential datatype
SOURCE
1365 subroutine effective_potential_setSupercell(eff_pot,comm,ncell,supercell) 1366 1367 !Arguments ------------------------------------ 1368 !scalars 1369 integer,intent(in) :: comm 1370 !arrays 1371 type(effective_potential_type),intent(inout) :: eff_pot 1372 integer,optional,intent(in) :: ncell(3) 1373 type(supercell_type),optional,intent(in) :: supercell 1374 !Local variables------------------------------- 1375 !scalar 1376 !arrays 1377 !character(len=500) :: msg 1378 1379 ! ************************************************************************* 1380 1381 !Checks 1382 if (.not.present(supercell).and..not.present(ncell)) then 1383 ABI_ERROR(' You should at least set ncell of supercell type') 1384 end if 1385 1386 call destroy_supercell(eff_pot%supercell) 1387 1388 if(present(supercell))then 1389 call copy_supercell(supercell,eff_pot%supercell) 1390 else 1391 call init_supercell(eff_pot%crystal%natom, (/ncell(1),0,0, 0,ncell(2),0, 0,0,ncell(3)/), & 1392 & eff_pot%crystal%rprimd,eff_pot%crystal%typat,eff_pot%crystal%xcart,& 1393 & eff_pot%crystal%znucl, eff_pot%supercell) 1394 end if 1395 1396 !Initialisation of new mpi over supercell 1397 call effective_potential_initmpi(eff_pot,comm) 1398 1399 end subroutine effective_potential_setSupercell
m_effective_potential/effective_potential_type [ Types ]
[ Top ] [ m_effective_potential ] [ Types ]
NAME
effective_potential_type
FUNCTION
datatype for a effective potential constructed.
SOURCE
104 type, public :: effective_potential_type 105 106 character(len=fnlen) :: name 107 ! Name of the molecule (CaTiO3,...) 108 109 type(crystal_t) :: crystal 110 ! crystal type 111 ! contains all information of the crystal 112 113 real(dp):: energy 114 ! Energy of the system (Hatree) 115 116 real(dp), allocatable :: fcart(:,:) 117 ! forces(3,natom) 118 ! initial cartesian forces of the system 119 120 real(dp) :: strten(6) 121 ! strten(6) 122 ! initial stresses (Ha/bohr^3) of the system 123 124 type(harmonics_terms_type) :: harmonics_terms 125 ! type with all information for harmonics terms 126 127 type(anharmonics_terms_type) :: anharmonics_terms 128 ! type with all information for anharmonics terms 129 130 type(polynomial_conf_type) :: confinement 131 ! type with all the information for the confinement 132 133 type(supercell_type) :: supercell 134 ! super cell type 135 ! Store all the information of the suppercell 136 137 logical :: has_anharmonicsTerms 138 ! True : the aharmonic part is present 139 140 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 141 ! This is for the parallelisation over the supercell 142 type(effpot_mpi_type) :: mpi_ifc 143 ! effpot_mpi_type with all the information for the IFC paralellisation 144 145 type(effpot_mpi_type) :: mpi_coeff 146 ! effpot_mpi_type with all the information for the polynomial coefficients paralellisation 147 148 end type effective_potential_type
m_effective_potential/effective_potential_writeAbiInput [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_writeAbiInput
FUNCTION
This routine print the effective potential into input of abinit We can also apply a strain to the structure
INPUTS
eff_pot<type(effective_potential_type)> = effective_potential datatype filename = the name of input file strain<(strain_type)> = optional,strain datatype if need to apply strain into rprim
OUTPUT
SOURCE
1998 subroutine effective_potential_writeAbiInput(eff_pot,filename,strain) 1999 2000 !Arguments ------------------------------------ 2001 !scalars 2002 type(strain_type),optional,intent(in) :: strain 2003 character(len=fnlen),optional,intent(in) :: filename 2004 !arrays 2005 type(effective_potential_type), intent(in) :: eff_pot 2006 !Local variables------------------------------- 2007 !scalar 2008 integer :: unit = 20 2009 character(len=500) :: msg 2010 character(len=fnlen) :: namefile 2011 !arrays 2012 real(dp) :: xred(3,eff_pot%crystal%natom) 2013 type(strain_type) :: strain_tmp 2014 2015 ! ************************************************************************ 2016 2017 if(present(strain)) then 2018 strain_tmp = strain 2019 else 2020 call strain_init(strain_tmp) 2021 end if 2022 2023 ! try to open the file 2024 if(present(filename)) then 2025 namefile=filename 2026 else 2027 if(eff_pot%name /='') then 2028 write(namefile,'(a)') 'structure_'//trim(eff_pot%name) 2029 else 2030 write(namefile,'(a)') 'structure' 2031 end if 2032 if (strain_tmp%name/='') then 2033 write(namefile,'(a,a,a,a,a,a,a)') trim(namefile)//"_"//trim(strain_tmp%name)//"_"//& 2034 & trim(itoa(strain_tmp%direction)),"_"//trim(ftoa(strain_tmp%delta)) 2035 end if 2036 2037 namefile=trim(namefile)//".in" 2038 2039 end if 2040 2041 call isfile(namefile,'new') 2042 2043 if (open_file(namefile,msg,unit=unit,form="formatted",status="new",action="write") /= 0) then 2044 ABI_ERROR(msg) 2045 end if 2046 2047 write(msg,'(a,a,a,a)')ch10,& 2048 & ' Generation of the input file in ',namefile,ch10 2049 call wrtout(ab_out,msg,'COLL') 2050 call wrtout(std_out,msg,'COLL') 2051 2052 write(unit,'("#Abinit Input for DFPT, this file contrains the keyword")') 2053 write(unit,'("#To run DFPT calculation of")') 2054 write(unit,'(a)') trim(eff_pot%name) 2055 if (strain_tmp%direction /= 0) then 2056 write(unit,'("# With a perturbation ")',advance="no") 2057 write(unit,'(a)',advance="no") trim(strain_tmp%name) 2058 write(unit,'(" in the direction : ")',advance="no") 2059 write(unit,'(a)',advance='no') trim(itoa(strain_tmp%direction)) 2060 write(unit,'(" with the deformation : ")',advance="no") 2061 write(unit,'(a)') trim(ftoa(strain_tmp%delta)) 2062 end if 2063 2064 write(unit,'("")') 2065 write(unit,'("ndtset 1 jdtset 1 2 3")') 2066 write(unit,'("")') 2067 write(unit,'("#DATASET1 GROUND STATE")') 2068 write(unit,'("tolwfr1 = 1d-15")') 2069 write(unit,'(" prtwf1 = 1")') 2070 write(unit,'(" nline1 = 5")') 2071 2072 write(unit,'("")') 2073 write(unit,'("#DATASET2 DDK PERTURBATION")') 2074 write(unit,'("getwfk2 = 1")') 2075 write(unit,'(" iscf2 = -3")') 2076 write(unit,'(" nline2 = 15")') 2077 write(unit,'("nnsclo2 = 5")') 2078 write(unit,'("kptopt2 = 2")') 2079 write(unit,'(" nqpt2 = 1")') 2080 write(unit,'(" qpt2 = 0 0 0 ")') 2081 write(unit,'("rfelfd2 = 2")') 2082 write(unit,'(" rfdir2 = 1 1 1 ")') 2083 write(unit,'("tolwfr2 = 1.0d-20 ")') 2084 write(unit,'(" prtwf2 = 1 ")') 2085 write(unit,'(" ")') 2086 2087 write(unit,'("#DATASET3 RF")') 2088 write(unit,'(" getddk3 = 2")') 2089 write(unit,'(" getwfk3 = 1")') 2090 write(unit,'(" iscf3 = 7")') 2091 write(unit,'(" kptopt3 = 2")') 2092 write(unit,'(" nqpt3 = 1")') 2093 write(unit,'(" qpt3 = 0 0 0")') 2094 write(unit,'(" rfphon3 = 1")') 2095 write(unit,'("rfatpol3 = 1 ")',advance='no') 2096 write(unit,'(a)') itoa(eff_pot%crystal%natom) 2097 write(unit,'(" rfelfd3 = 3")') 2098 write(unit,'(" rfstrs3 = 3")') 2099 write(unit,'(" rfdir3 = 1 1 1")') 2100 write(unit,'(" tolvrs3 = 1.0d-8")') 2101 write(unit,'("")') 2102 2103 write(unit,'("#STRUCTURE")') 2104 write(unit,'(" natom = ")',advance='no') 2105 write(unit,'(a)') itoa(eff_pot%crystal%natom) 2106 write(unit,'(" znucl =")',advance='no') 2107 write(unit,'(10(F4.0))') (eff_pot%crystal%znucl) 2108 write(unit,'("ntypat = ")',advance='no') 2109 write(unit,'(a)') itoa(eff_pot%crystal%ntypat) 2110 write(unit,'(" typat = ")',advance='no') 2111 write(unit,'(10(I2))') (eff_pot%crystal%typat) 2112 write(unit,'(" acell = 1 1 1")') 2113 write(unit,'(" rprim ")') 2114 write(unit,'(3(F20.10))') (matmul(eff_pot%crystal%rprimd,strain%strain)) 2115 write(unit,'(" xred ")') 2116 call xcart2xred(eff_pot%crystal%natom,eff_pot%crystal%rprimd,eff_pot%crystal%xcart,xred) 2117 write(unit,'(3(F15.10))') (xred) 2118 write(unit,'(" ")') 2119 2120 write(unit,'("#SCF")') 2121 write(unit,'(" ecut = ")') 2122 write(unit,'("pawecutdg = ")') 2123 write(unit,'(" ecutsm = ")') 2124 write(unit,'(" tolvrs = ")') 2125 write(unit,'(" nband = ")') 2126 write(unit,'(" ixc = ")') 2127 write(unit,'(" occopt = ")') 2128 write(unit,'(" nstep = ")') 2129 write(unit,'(" kptopt = ")') 2130 write(unit,'(" ngkpt = ")') 2131 write(unit,'("")') 2132 write(unit,'(" prtwf 0 prtden 0 prtdos 0")') 2133 2134 close(unit) 2135 2136 end subroutine effective_potential_writeAbiInput
m_effective_potential/effective_potential_writeNETCDF [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_writeNETCDF
FUNCTION
This routine print the effective potential into netcdf format Several options are available
INPUTS
filename = the name of output file
eff_pot = datatype contains the effective potential
option = option for the format of the xml file
1 print the xml for a system
OUTPUT
SOURCE
3598 subroutine effective_potential_writeNETCDF(eff_pot,option,filename) 3599 3600 !Arguments ------------------------------------ 3601 !scalars 3602 integer, intent(in) :: option 3603 character(len=fnlen),optional,intent(in) :: filename 3604 !arrays 3605 type(effective_potential_type), intent(in) :: eff_pot 3606 3607 !Local variables------------------------------- 3608 !scalar 3609 integer :: amu_id,bec_id,ifccell_id,epsinf_id,elastic_id 3610 integer :: ifc_id,ifcs_id,natom_id,ntypat_id,nrpt_id,npsp_id,typat_id 3611 integer :: six_id,two_id,xyz_id,znucl_id 3612 integer :: ncerr,ncid,npsp 3613 integer :: dimCids(2),dimEids(2),dimIids(6),dimPids(1),dimRids(2),dimXids(2) 3614 integer :: etotal_id,rprimd_id,xcart_id 3615 character(len=500) :: msg 3616 character(len=fnlen) :: namefile 3617 !arrays 3618 real(dp) :: strain(9,6) 3619 3620 ! ************************************************************************* 3621 3622 #if defined HAVE_NETCDF 3623 3624 strain(:,1) = (/1,0,0,0,0,0,0,0,0/) 3625 strain(:,2) = (/0,0,0,0,1,0,0,0,0/) 3626 strain(:,3) = (/0,0,0,0,0,0,0,0,1/) 3627 strain(:,4) = half*(/0,0,0,0,0,1,0,1,0/) 3628 strain(:,5) = half*(/0,0,1,0,0,0,1,0,0/) 3629 strain(:,6) = half*(/0,1,0,1,0,0,0,0,0/) 3630 3631 !Print only the reference system in xml format 3632 if (option == 1) then 3633 3634 if(present(filename)) then 3635 namefile=filename 3636 else 3637 namefile='ref.nc' 3638 end if 3639 3640 call isfile(namefile,'new') 3641 3642 write(msg,'(a,a,a)')ch10,& 3643 & ' Generation of the xml file for the reference structure in ',namefile 3644 3645 call wrtout(ab_out,msg,'COLL') 3646 call wrtout(std_out,msg,'COLL') 3647 3648 ! 1. Create netCDF file 3649 ncerr = nf90_create(path=trim(namefile),cmode=NF90_CLOBBER, ncid=ncid) 3650 NCF_CHECK_MSG(ncerr,"create netcdf history file") 3651 3652 ! 2. Define dimensions 3653 ncerr = nf90_def_dim(ncid,"natom",eff_pot%crystal%natom,natom_id) 3654 NCF_CHECK_MSG(ncerr," define dimension natom") 3655 3656 ncerr = nf90_def_dim(ncid,"ntypat",eff_pot%crystal%ntypat,ntypat_id) 3657 NCF_CHECK_MSG(ncerr," define dimension ntypat") 3658 3659 ncerr = nf90_def_dim(ncid,"nrpt",eff_pot%harmonics_terms%ifcs%nrpt,nrpt_id) 3660 NCF_CHECK_MSG(ncerr," define dimension ntypat") 3661 3662 ncerr = nf90_def_var(ncid, "typat", NF90_INT, natom_id, typat_id) 3663 NCF_CHECK_MSG(ncerr," define variable typat") 3664 3665 npsp = size(eff_pot%crystal%znucl) 3666 if (npsp /= eff_pot%crystal%ntypat) then 3667 ABI_WARNING("HIST file does not support alchemical mixing!") 3668 end if 3669 ncerr = nf90_def_dim(ncid,"npsp",npsp,npsp_id) 3670 NCF_CHECK_MSG(ncerr," define dimension npsp") 3671 3672 ncerr = nf90_def_var(ncid, "znucl", NF90_DOUBLE, npsp_id, znucl_id) 3673 NCF_CHECK_MSG(ncerr," define variable znucl") 3674 3675 ncerr = nf90_def_dim(ncid,"xyz",3,xyz_id) 3676 NCF_CHECK_MSG(ncerr," define dimension xyz") 3677 3678 ncerr = nf90_def_dim(ncid,"six",6,six_id) 3679 NCF_CHECK_MSG(ncerr," define dimension six") 3680 3681 ncerr = nf90_def_dim(ncid,"two",2,two_id) 3682 NCF_CHECK_MSG(ncerr," define dimension two") 3683 3684 ! Dimensions for xcart,xred,fcart,gred and vel 3685 dimXids = (/ xyz_id, natom_id /) 3686 ! Dimensions for rprimd 3687 dimRids = (/ xyz_id, xyz_id /) 3688 ! Dimensions for ifc 3689 dimIids = (/ 2, xyz_id, natom_id, xyz_id, natom_id, nrpt_id /) 3690 ! Dimensions for position 3691 dimPids = (/ xyz_id /) 3692 ! Dimension for elastic constant 3693 dimEids = (/six_id,six_id/) 3694 ! Dimension for cell 3695 dimCids = (/nrpt_id,3/) 3696 3697 ! 3. Define variables and their attributes (units and mnemonics) 3698 call ab_define_var(ncid, (/1/), etotal_id, NF90_DOUBLE,& 3699 & "energy","Energy of the reference structure","Ha" ) 3700 3701 call ab_define_var(ncid, dimRids, rprimd_id, NF90_DOUBLE,& 3702 & "rprimd","Real space PRIMitive translations, Dimensional","bohr" ) 3703 3704 call ab_define_var(ncid, dimRids, epsinf_id, NF90_DOUBLE,& 3705 & "epsilon_inf","Dielectric tensor, Dimensional","epsilon_inf" ) 3706 3707 call ab_define_var(ncid, dimEids, elastic_id, NF90_DOUBLE,& 3708 & "elastic","Elastic Constants, Dimensional","Ha" ) 3709 3710 call ab_define_var(ncid, dimRids, bec_id, NF90_DOUBLE,& 3711 & "bec","Born Effective Charges, Dimensional","abs(e)" ) 3712 3713 call ab_define_var(ncid, dimXids, xcart_id, NF90_DOUBLE,& 3714 & "xcart","vectors (X) of atom positions in CARTesian coordinates","bohr" ) 3715 3716 call ab_define_var(ncid, dimIids, ifcs_id, NF90_DOUBLE,& 3717 & "IFCs","Interatomic Forces Constantes in real spaces (short range), Dimensional","Hatree/bohr**2" ) 3718 3719 call ab_define_var(ncid, dimIids, ifc_id, NF90_DOUBLE,& 3720 & "IFC","Interatomic Forces Constantes in real spaces (total range), Dimensional","Hatree/bohr**2" ) 3721 3722 call ab_define_var(ncid, dimCids, ifccell_id, NF90_DOUBLE,& 3723 & "cell","cell for the ifc, Dimensional","Dimensionless" ) 3724 3725 call ab_define_var(ncid, [ntypat_id], amu_id, NF90_DOUBLE,& 3726 & "amu","Masses of each type of atom in atomic mass units", "" ) 3727 3728 ! 4. End define mode 3729 ncerr = nf90_enddef(ncid) 3730 NCF_CHECK_MSG(ncerr," end define mode") 3731 3732 ! 5. Write variables 3733 ncerr = nf90_put_var(ncid,etotal_id, eff_pot%energy) 3734 NCF_CHECK_MSG(ncerr," write variable energy") 3735 3736 ncerr = nf90_put_var(ncid,rprimd_id, eff_pot%crystal%rprimd) 3737 NCF_CHECK_MSG(ncerr," write variable rprimd") 3738 3739 ncerr = nf90_put_var(ncid,epsinf_id, eff_pot%harmonics_terms%epsilon_inf) 3740 NCF_CHECK_MSG(ncerr," write variable epsilon_inf") 3741 3742 ncerr = nf90_put_var(ncid,elastic_id , eff_pot%harmonics_terms%elastic_constants) 3743 NCF_CHECK_MSG(ncerr," write variable elastic_constant") 3744 3745 ncerr = nf90_put_var(ncid, bec_id, eff_pot%harmonics_terms%zeff) 3746 NCF_CHECK_MSG(ncerr," write variable bec") 3747 3748 ncerr = nf90_put_var(ncid,xcart_id, eff_pot%crystal%xcart) 3749 NCF_CHECK_MSG(ncerr," write variable xcart") 3750 3751 ncerr = nf90_put_var(ncid,ifccell_id, eff_pot%harmonics_terms%ifcs%cell) 3752 NCF_CHECK_MSG(ncerr," write variable cell") 3753 3754 ncerr = nf90_put_var(ncid,ifcs_id, eff_pot%harmonics_terms%ifcs%short_atmfrc) 3755 NCF_CHECK_MSG(ncerr," write variable short ifc") 3756 3757 ncerr = nf90_put_var(ncid,ifc_id, eff_pot%harmonics_terms%ifcs%atmfrc) 3758 NCF_CHECK_MSG(ncerr," write variable total ifc") 3759 3760 3761 ! 6. Close NetCDF file 3762 ncerr = nf90_close(ncid) 3763 NCF_CHECK_MSG(ncerr," close netcdf history file") 3764 end if 3765 3766 #endif 3767 end subroutine effective_potential_writeNETCDF
m_effective_potential/effective_potential_writeXML [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
effective_potential_writeXML
FUNCTION
This routine print the effective potential into xml format Several options are available
INPUTS
filename = the name of output file eff_pot<type(effective_potential_type)> = effective_potential datatype option = 0 Do nothing = 1 Generate the XML file with: - The system definition and the model (Harmonic + Anharmonic) = 2 Generate two XML files with: - The system definition and the model (Harmonic) - The model (Anharmonic) = 3 Generate one XML files with: - The system definition and the model (Harmonic) = 4 Generate one XML files with: - The model (Anharmonic)
OUTPUT
SOURCE
1679 subroutine effective_potential_writeXML(eff_pot,option,filename,prt_dipdip) 1680 1681 !Arguments ------------------------------------ 1682 !scalars 1683 integer, intent(in) :: option 1684 character(len=fnlen),optional,intent(in) :: filename 1685 logical,optional,intent(in) :: prt_dipdip 1686 !arrays 1687 type(effective_potential_type), intent(in) :: eff_pot 1688 1689 !Local variables------------------------------- 1690 !scalar 1691 integer :: ii,ia,ib,jj 1692 integer :: iqpt,irpt,mu,nu 1693 integer :: unit_xml 1694 character(len=500) :: msg 1695 character(len=fnlen) :: namefile 1696 character(len=10) :: natom 1697 logical :: new_file,need_prtdipdip 1698 !arrays 1699 real(dp) :: strain(9,6) 1700 1701 ! ************************************************************************* 1702 1703 strain(:,1) = (/1,0,0,0,0,0,0,0,0/) 1704 strain(:,2) = (/0,0,0,0,1,0,0,0,0/) 1705 strain(:,3) = (/0,0,0,0,0,0,0,0,1/) 1706 strain(:,4) = half*(/0,0,0,0,0,1,0,1,0/) 1707 strain(:,5) = half*(/0,0,1,0,0,0,1,0,0/) 1708 strain(:,6) = half*(/0,1,0,1,0,0,0,0,0/) 1709 1710 unit_xml = get_unit() 1711 need_prtdipdip = .TRUE. 1712 if(present(prt_dipdip)) need_prtdipdip = prt_dipdip 1713 !Print only the reference system in xml format 1714 if (option == 1 .or. option == 2 .or. option ==3) then 1715 1716 ! convert natom in character 1717 write (natom,'(I9)') eff_pot%crystal%natom 1718 1719 ! Compute the name of the XML file 1720 if(present(filename)) then 1721 namefile=replace(trim(filename),".out","") 1722 select case(option) 1723 case(1) 1724 namefile=trim(namefile)//"_model.xml" 1725 case(2) 1726 namefile=trim(namefile)//"_sys.xml" 1727 case(3) 1728 namefile=trim(namefile)//"_sys.xml" 1729 end select 1730 else 1731 namefile='system.xml' 1732 end if 1733 1734 call isfile(namefile,'new') 1735 1736 if (open_file(namefile,msg,unit=unit_xml,form="formatted",& 1737 & status="new",action="write") /= 0) then 1738 ABI_ERROR(msg) 1739 end if 1740 1741 write(msg,'(a,a,a)')ch10,& 1742 & ' Generation of the xml file for the model in ',trim(namefile) 1743 1744 call wrtout(ab_out,msg,'COLL') 1745 call wrtout(std_out,msg,'COLL') 1746 1747 ! Write header 1748 WRITE(unit_xml,'("<?xml version=""1.0"" ?>")') 1749 WRITE(unit_xml,'("<System_definition>")') 1750 1751 WRITE(unit_xml,'(" <energy>")') 1752 WRITE(unit_xml,'(E23.14)') (eff_pot%energy) 1753 WRITE(unit_xml,'(" </energy>")') 1754 1755 WRITE(unit_xml,'(" <unit_cell units=""bohrradius"">")') 1756 do mu=1,3 1757 WRITE(unit_xml,'(3(E23.14))') (eff_pot%crystal%rprimd(mu,nu),nu=1,3) 1758 end do 1759 WRITE(unit_xml,'(" </unit_cell>")') 1760 1761 WRITE(unit_xml,'(" <epsilon_inf units=""epsilon0"">")') 1762 WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%epsilon_inf) 1763 WRITE(unit_xml,'(" </epsilon_inf>")') 1764 1765 WRITE(unit_xml,'(" <elastic units=""hartree"">")') 1766 WRITE(unit_xml,'(6(E23.14))') (eff_pot%harmonics_terms%elastic_constants) 1767 WRITE(unit_xml,'(" </elastic>")') 1768 1769 do ia=1,eff_pot%crystal%natom 1770 WRITE(unit_xml,'(" <atom mass=""",1F10.5,""" massunits=""atomicmassunit"">")') & 1771 & eff_pot%crystal%amu(eff_pot%crystal%typat(ia)) 1772 WRITE(unit_xml,'(" <position units=""bohrradius"">")') 1773 WRITE(unit_xml,'(3(E23.14))') (eff_pot%crystal%xcart(:,ia)) 1774 WRITE(unit_xml,'(" </position>")') 1775 WRITE(unit_xml,'(" <borncharge units=""abs(e)"">")') 1776 WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%zeff(:,:,ia)) 1777 WRITE(unit_xml,'(" </borncharge>")') 1778 WRITE(unit_xml,'(" </atom>")') 1779 end do 1780 ! Print the Ifc short range for each cell data is array 3*natom*3*natom 1781 ! [ [x1 x2 ....] 1782 ! [y1 y2 ....] for atom 1 1783 ! [z1 z2 ....] 1784 ! [x1 x2 ....] 1785 ! [y1 y2 ....] for atom 2 1786 ! [z1 z2 ....] 1787 ! .... ] 1788 ! Warning : The IFC are print in other order 1,mu,ia,nu,ib which is not fortran way 1789 ! We do like that to because the previous script was in python 1790 ! When you read ifc from XML file with fotran, you have to tranpose the matrix 1791 ! 1792 do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt 1793 if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol20)) then 1794 WRITE(unit_xml,'(" <local_force_constant units=""hartree/bohrradius**2"">")') 1795 WRITE(unit_xml,'(" <data>")') 1796 do ia=1,eff_pot%crystal%natom 1797 do mu=1,3 1798 do ib=1,eff_pot%crystal%natom 1799 do nu=1,3 1800 WRITE(unit_xml,'(e22.14)', advance="no")& 1801 & (eff_pot%harmonics_terms%ifcs%short_atmfrc(mu,ia,nu,ib,irpt)) 1802 end do 1803 end do 1804 WRITE(unit_xml,'(a)')'' 1805 end do 1806 end do 1807 WRITE(unit_xml,'(" </data>")') 1808 WRITE(unit_xml,'(" <cell>")') 1809 WRITE(unit_xml,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt)) 1810 WRITE(unit_xml,'(" </cell>")') 1811 WRITE(unit_xml,'(" </local_force_constant>")') 1812 end if 1813 ! Print the IFC total for each cell, data is array 3*natom*3*natom 1814 ! [ [x1 x2 ....] 1815 ! [y1 y2 ....] for atom 1 1816 ! [z1 z2 ....] 1817 ! [x1 x2 ....] 1818 ! [y1 y2 ....] for atom 2 1819 ! [z1 z2 ....] 1820 ! .... ] 1821 if(need_prtdipdip)then 1822 if(all(abs(eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt))<tol20)) then 1823 if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol20)) then 1824 write(msg, '(a,a,a,a)' )& 1825 & ' There is no total range but short range in your effective potential',ch10,& 1826 & 'Action: contact abinit group' 1827 ABI_BUG(msg) 1828 end if 1829 else 1830 WRITE(unit_xml,'(" <total_force_constant units=""hartree/bohrradius**2"">")') 1831 WRITE(unit_xml,'(" <data>")') 1832 do ia=1,eff_pot%crystal%natom 1833 do mu=1,3 1834 do ib=1,eff_pot%crystal%natom 1835 do nu=1,3 1836 WRITE(unit_xml,'(e22.14)', advance="no")& 1837 & (eff_pot%harmonics_terms%ifcs%atmfrc(mu,ia,nu,ib,irpt)) 1838 end do 1839 end do 1840 WRITE(unit_xml,'(a)')'' 1841 end do 1842 end do 1843 WRITE(unit_xml,'(" </data>")') 1844 WRITE(unit_xml,'(" <cell>")') 1845 WRITE(unit_xml,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt)) 1846 WRITE(unit_xml,'(" </cell>")') 1847 WRITE(unit_xml,'(" </total_force_constant>")') 1848 end if 1849 end if 1850 end do 1851 1852 do iqpt=1,eff_pot%harmonics_terms%nqpt 1853 WRITE(unit_xml,'(" <phonon>")') 1854 WRITE(unit_xml,'(" <qpoint units=""2pi*G0"">")') 1855 WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%qpoints(:,iqpt)) 1856 WRITE(unit_xml,'(" </qpoint>")') 1857 WRITE(unit_xml,'(" <frequencies units=""reciprocal cm"">")') 1858 WRITE(unit_xml,'(3(e22.14))') (eff_pot%harmonics_terms%phfrq(:,iqpt)) 1859 WRITE(unit_xml,'(" </frequencies>")') 1860 WRITE(unit_xml,'(" <dynamical_matrix units=""hartree/bohrradius**2"">")') 1861 do ia=1,eff_pot%crystal%natom 1862 do mu=1,3 1863 do ib = 1,eff_pot%crystal%natom 1864 do nu=1,3 1865 WRITE(unit_xml,'(e22.14)',advance='no')(eff_pot%harmonics_terms%dynmat(1,nu,ib,mu,ia,iqpt)) 1866 end do 1867 end do 1868 WRITE(unit_xml,'(a)')'' 1869 end do 1870 end do 1871 WRITE(unit_xml,'(" </dynamical_matrix>")') 1872 WRITE(unit_xml,'(" </phonon>")') 1873 end do 1874 ! if phonon/forces strain is computed 1875 jj = 1 1876 do ii = 1,6 1877 WRITE(unit_xml,'(" <strain_coupling voigt=""",I2,""">")') ii-1 1878 WRITE(unit_xml,'(" <strain>")') 1879 WRITE(unit_xml,'(6(e12.4))') (strain(:,jj)) 1880 WRITE(unit_xml,'(" </strain>")') 1881 WRITE(unit_xml,'(" <correction_force units=""hartree/bohrradius"">")') 1882 do ia=1,eff_pot%crystal%natom 1883 do mu=1,3 1884 WRITE(unit_xml,'(e22.14)', advance="no")& 1885 & (eff_pot%harmonics_terms%strain_coupling(ii,mu,ia)) 1886 end do 1887 WRITE(unit_xml,'(a)')'' 1888 end do 1889 WRITE(unit_xml,'(" </correction_force>")') 1890 if (eff_pot%has_anharmonicsTerms)then 1891 if (eff_pot%anharmonics_terms%has_elastic3rd) then 1892 WRITE(unit_xml,'(" <elastic3rd units=""hartree"">")') 1893 WRITE(unit_xml,'(6(E23.14))') (eff_pot%anharmonics_terms%elastic3rd(ii,:,:)) 1894 WRITE(unit_xml,'(" </elastic3rd>")') 1895 end if 1896 if (eff_pot%anharmonics_terms%has_elastic_displ) then 1897 WRITE(unit_xml,'(" <correction_strain_force units=""hartree/bohrradius"">")') 1898 do ia=1,eff_pot%crystal%natom 1899 do mu=1,3 1900 do nu=1,6 1901 WRITE(unit_xml,'(e22.14)', advance="no")& 1902 & (eff_pot%anharmonics_terms%elastic_displacement(ii,nu,mu,ia)) 1903 end do 1904 end do 1905 WRITE(unit_xml,'(a)')'' 1906 end do 1907 WRITE(unit_xml,'(" </correction_strain_force>")') 1908 end if 1909 if (eff_pot%anharmonics_terms%has_strain_coupling) then 1910 do irpt=1,eff_pot%anharmonics_terms%phonon_strain(ii)%nrpt 1911 WRITE(unit_xml,'(" <correction_force_constant units=""hartree/bohrradius**2"">")') 1912 WRITE(unit_xml,'(" <data>")') 1913 do ia=1,eff_pot%crystal%natom 1914 do mu=1,3 1915 do ib=1,eff_pot%crystal%natom 1916 do nu=1,3 1917 WRITE(unit_xml,'(e22.14)', advance="no")& 1918 & (eff_pot%anharmonics_terms%phonon_strain(ii)%atmfrc(mu,ia,nu,ib,irpt)) 1919 end do 1920 end do 1921 WRITE(unit_xml,'(a)')'' 1922 end do 1923 end do 1924 WRITE(unit_xml,'(" </data>")') 1925 WRITE(unit_xml,'(" <cell>")') 1926 WRITE(unit_xml,'(3(I4))') (eff_pot%anharmonics_terms%phonon_strain(ii)%cell(:,irpt)) 1927 WRITE(unit_xml,'(" </cell>")') 1928 WRITE(unit_xml,'(" </correction_force_constant>")') 1929 end do 1930 end if!End if has_straincouplitn 1931 end if!end Hasstrain_coupling 1932 WRITE(unit_xml,'(" </strain_coupling>")') 1933 jj = jj + 1 1934 end do!end mu 1935 1936 if(option /=1) WRITE(unit_xml,'("</System_definition>")') 1937 ! Close file 1938 CLOSE(unit_xml) 1939 1940 end if!end option 1941 1942 !Print the coefficients into XML file 1943 new_file = .FALSE. 1944 if (option==1 .or. option == 2 .or. option==4) then 1945 ! Compute the name of the XML file 1946 if(present(filename)) then 1947 namefile=replace(trim(filename),".out","") 1948 select case(option) 1949 case(1) 1950 new_file = .FALSE. 1951 namefile=trim(namefile)//"_model.xml" 1952 case(2) 1953 new_file = .TRUE. 1954 namefile=trim(namefile)//"_coeffs.xml" 1955 case(4) 1956 new_file = .TRUE. 1957 namefile=trim(namefile)//"_coeffs.xml" 1958 end select 1959 else 1960 namefile='coeffs.xml' 1961 end if 1962 1963 if(eff_pot%anharmonics_terms%ncoeff > 0) then 1964 call polynomial_coeff_writeXML(eff_pot%anharmonics_terms%coefficients,& 1965 & eff_pot%anharmonics_terms%ncoeff,namefile,unit=unit_xml,& 1966 & newfile=new_file) 1967 end if 1968 end if!end option 1969 1970 if(option==1)then 1971 ! add the end of the file in the case option==1 1972 open(unit=unit_xml,file=namefile,position="append") 1973 WRITE(unit_xml,'("</System_definition>")') 1974 ! Close file 1975 CLOSE(unit_xml) 1976 end if 1977 1978 end subroutine effective_potential_writeXML
m_effective_potential/equal [ Functions ]
[ Top ] [ m_effective_potential ] [ Functions ]
NAME
equal
FUNCTION
compare two effective potential
INPUTS
e1<type(effective_potential_type)> = effective_potential datatype e2<type(effective_potential_type)> = effective_potential datatype
OUTPUT
SOURCE
3007 pure function effective_potential_compare(e1,e2) result (res) 3008 3009 !Arguments ------------------------------------ 3010 type(effective_potential_type), intent(in) :: e1,e2 3011 logical :: res 3012 ! ************************************************************************* 3013 res = .false. 3014 if(e1%crystal%natom==e2%crystal%natom.and.& 3015 & e1%harmonics_terms%ifcs%nrpt==e2%harmonics_terms%ifcs%nrpt.and.& 3016 & e1%crystal%ntypat==e2%crystal%ntypat.and.& 3017 & e1%harmonics_terms%nqpt==e2%harmonics_terms%nqpt.and.& 3018 & abs(e1%energy-e2%energy)<tol16.and.& 3019 & abs(e1%crystal%ucvol-e2%crystal%ucvol)<tol16) then 3020 res = .true. 3021 end if 3022 3023 end function effective_potential_compare