TABLE OF CONTENTS


ABINIT/m_effective_potential [ Modules ]

[ Top ] [ 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-2018 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

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_effective_potential
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_errors
33  use m_abicore
34  use m_strain
35  use m_ifc
36  use m_supercell
37  use m_phonons
38  use m_ddb
39  use m_polynomial_conf
40  use m_polynomial_coeff
41  use m_anharmonics_terms
42  use m_harmonics_terms
43  use m_xmpi
44  use m_ewald
45  use m_nctk
46 #if defined HAVE_NETCDF
47  use netcdf
48 #endif
49 
50  use m_fstrings,       only : replace, ftoa, itoa, int2char4
51  use m_io_tools,       only : open_file,get_unit
52  use m_dtfil,          only : isfile
53  use m_symtk,          only : matr3inv
54  use m_effpot_mpi,     only : effpot_mpi_init,effpot_mpi_type,effpot_mpi_free
55  use m_abihist,        only : abihist
56  use m_special_funcs,  only : factorial
57  use m_copy,           only : alloc_copy
58  use m_geometry,       only : fred2fcart,fcart2fred, xcart2xred, xred2xcart, metric
59  use m_crystal,        only : crystal_t, crystal_init, crystal_free,crystal_print
60  use m_anaddb_dataset, only : anaddb_dataset_type, anaddb_dtset_free, outvars_anaddb, invars9
61  use m_multibinit_dataset, only : multibinit_dtset_type
62  use m_dynmat,         only : make_bigbox,q0dy3_apply, q0dy3_calc, dfpt_phfrq
63 
64  implicit none
65 
66  public :: effective_potential_distributeResidualForces
67  public :: effective_potential_evaluate
68  public :: effective_potential_free
69  public :: effective_potential_freeCoeffs
70  public :: effective_potential_freempi
71  public :: effective_potential_generateDipDip
72  public :: effective_potential_getDisp
73  public :: effective_potential_init
74  public :: effective_potential_initmpi
75  public :: effective_potential_print
76  public :: effective_potential_printSupercell
77  public :: effective_potential_setCoeffs
78  public :: effective_potential_setConfinement
79  public :: effective_potential_setElastic3rd
80  public :: effective_potential_setElastic4th
81  public :: effective_potential_setElasticDispCoupling
82  public :: effective_potential_setStrainPhononCoupling
83  public :: effective_potential_setSupercell
84  public :: effective_potential_writeAbiInput
85  public :: effective_potential_writeXML
86  !AM_EXPERIMENTAL
87  public :: effective_potential_computeGradient
88 ! public :: effective_potential_effpot2ddb
89  ! public :: effective_potential_printPDOS
90  public :: effective_potential_checkDEV
91  public :: effective_potential_writeNETCDF
92  public :: OPERATOR(==)
93  !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

PARENTS

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

3519  subroutine effective_potential_checkDEV(eff_pot,hist,natom,ntime)
3520 
3521 
3522 !This section has been created automatically by the script Abilint (TD).
3523 !Do not modify the following lines by hand.
3524 #undef ABI_FUNC
3525 #define ABI_FUNC 'effective_potential_checkDEV'
3526 !End of the abilint section
3527 
3528  implicit none
3529 
3530 !Arguments ------------------------------------
3531 !scalars
3532  integer, intent(in) :: natom,ntime
3533 !arrays
3534  type(effective_potential_type),intent(in) :: eff_pot
3535  type(abihist),intent(in) :: hist
3536 !Local variables-------------------------------
3537 !scalar
3538  integer :: ii,jj,ia,mu,npt,istep
3539 ! integer :: ifirst
3540  real(dp):: energy,delt,delta,ucvol
3541  !arrays
3542  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),mat_def(3,3),identity(3,3)
3543  real(dp):: fcart(3,natom),fred(3,natom),strten(6),rprimd(3,3)
3544  real(dp):: rprimd_def(3,3),rprimd_ref(3,3),deltalist(5)
3545  real(dp):: disp(3,natom),disp_red(3,natom),strain(6),du_delta(6,3,natom),diff(5)
3546  real(dp),allocatable :: xred(:,:)
3547  integer,parameter :: alpha(9)=(/1,2,3,3,3,2,2,1,1/),beta(9)=(/1,2,3,2,1,1,3,3,2/)
3548  character(len=500) :: msg
3549 
3550 ! *************************************************************************
3551 
3552  !Do some checks
3553  if(ntime /= hist%mxhist)then
3554    write(msg,'(a)')'ntime is not correct'
3555    MSG_BUG(msg)
3556  end if
3557 
3558  if(natom /= size(hist%xred,2)) then
3559    write(msg,'(a)')'natom is not correct'
3560    MSG_BUG(msg)
3561  end if
3562 
3563 
3564  ABI_ALLOCATE(xred,(3,natom))
3565  xred = zero
3566 
3567 !option 1 => set the reference for the test
3568 ! call xcart2xred(eff_pot%supercell%natom,eff_pot%supercell%rprimd,&
3569 !&                eff_pot%supercell%xcart,xred)
3570 ! rprimd =  eff_pot%supercell%rprimd
3571 
3572 !option 2 => set a specific step for the test
3573  istep = 4
3574  xred = hist%xred(:,:,istep)
3575  rprimd =  hist%rprimd(:,:,istep)
3576 
3577  rprimd_ref =  eff_pot%supercell%rprimd
3578  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
3579 
3580  npt=5
3581  delta = 0.001
3582  deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/)
3583  strain = zero
3584 
3585    do ia=1,natom
3586      do mu=1,3
3587        write(std_out,*) "atm: ",ia," dir: ",mu
3588        do ii=1,npt
3589          delt = deltalist(ii)
3590 
3591 !        Get the initial displacement
3592          call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3593 &                                         eff_pot%supercell%rprimd,1,xred_hist=xred,&
3594 &                                         xcart_ref=eff_pot%supercell%xcart,&
3595 &                                         compute_displacement = .true.,compute_duDelta = .true.)
3596 
3597 !        Add the delta
3598          call xcart2xred(natom, rprimd, disp, disp_red)
3599          disp_red(mu,ia) = disp_red(mu,ia) + delt
3600          call xred2xcart(natom, rprimd, disp, disp_red)
3601 
3602          call effective_potential_evaluate(eff_pot,energy,fcart,fred,strten,natom,rprimd,&
3603 &                                          xred=xred,du_delta=du_delta,&
3604 &                                          displacement=disp,compute_anharmonic=.true.,verbose=.false.)
3605          diff(ii) = energy
3606 
3607        end do
3608 
3609 !  Get the initial displacement
3610    call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3611 &                                   eff_pot%supercell%rprimd,1,xred_hist=xred,&
3612 &                                   xcart_ref=eff_pot%supercell%xcart,&
3613 &                                   compute_displacement = .true.,compute_duDelta = .true.)
3614 
3615    call effective_potential_evaluate(eff_pot,energy,fcart,fred,strten,natom,rprimd,&
3616 &                                    xred=xred,du_delta=du_delta,&
3617 &                                    displacement=disp,compute_anharmonic=.true.,verbose=.false.)
3618 
3619    write(std_out,*) "Analyti:",fred(mu,ia)
3620    write(std_out,*) "FD     :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta)
3621    write(std_out,*) "Diff(%):",abs(100*(fred(mu,ia)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))&
3622 &             / (12*delta) )) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) ))
3623 
3624  end do
3625 end do
3626 
3627 
3628 !  Fill the identity matrix
3629 identity = zero
3630 forall(ii=1:3)identity(ii,ii)=1
3631 
3632  npt=5
3633  delta = 0.0005
3634  deltalist = (/-2*delta,-delta,real(0.0,dp),delta,2*delta/)
3635 
3636  do jj=1,6
3637    write(std_out,*) "strain ",jj
3638    do ii=1,npt
3639      strain = zero
3640      delt = deltalist(ii)
3641      mat_def = zero
3642      strain(jj) = strain(jj) + delt
3643 
3644      mat_def(alpha(jj),beta(jj)) = mat_def(alpha(jj),beta(jj)) + half * strain(jj)
3645      mat_def(beta(jj),alpha(jj)) = mat_def(beta(jj),alpha(jj)) + half * strain(jj)
3646 
3647      rprimd_def =  matmul(identity(:,:)+mat_def(:,:),rprimd)
3648 
3649 ! The two options should give the same result
3650 ! Option 1 => compute the disps and provide them to evaluate
3651 !      call effective_potential_getDisp(disp,du_delta,natom,rprimd_def,&
3652 ! &                                     rprimd_ref,1,xred_hist=xred,&
3653 ! &                                     xcart_ref=eff_pot%supercell%xcart,&
3654 ! &                                     compute_displacement = .true.,compute_duDelta = .true.)
3655 
3656 !      call effective_potential_evaluate(eff_pot,energy,fcart,fred,strten,natom,rprimd_def,&
3657 ! &                                      xred=xred,du_delta=du_delta,&
3658 ! &                                      displacement=disp,strain=strain,&
3659 ! &                                      compute_anharmonic=.true.,verbose=.false.)
3660 
3661 !   Option 2 => compute the disps within evaluate
3662     call effective_potential_evaluate(eff_pot,energy,fcart,fred,strten,natom,rprimd_def,&
3663 &                                     xred=xred,compute_anharmonic=.true.,verbose=.false.)
3664 
3665 
3666 
3667      diff(ii) = energy
3668 
3669    end do
3670 
3671 !  The two options should give the same result
3672 !  Option 1 => compute the disps and provide them to evaluate
3673 !    call effective_potential_getDisp(disp,du_delta,natom,rprimd,&
3674 ! &                                   rprimd_ref,1,xred_hist=xred,&
3675 ! &                                   xcart_ref=eff_pot%supercell%xcart,&
3676 ! &                                   compute_displacement = .true.,compute_duDelta = .true.)
3677 
3678 !    call effective_potential_evaluate(eff_pot,energy,fcart,fred,strten,natom,rprimd,&
3679 ! &                                    xred=xred,du_delta=du_delta,&
3680 ! &                                    displacement=disp,&
3681 ! &                                    compute_anharmonic=.true.,verbose=.false.)
3682 
3683 !  Option 2 => compute the disps within evaluate
3684    call effective_potential_evaluate(eff_pot,energy,fcart,fred,strten,natom,rprimd,&
3685 &                                    xred=xred,compute_anharmonic=.true.,verbose=.false.)
3686 
3687  write(std_out,*) "Analyti:",strten(jj)
3688  write(std_out,*) "FD     :",(-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol
3689  write(std_out,*) "Diff(%):",abs(100*(strten(jj)-((-diff(5)+8*diff(4)-8*diff(2)+diff(1))&
3690 &   / (12*delta) / ucvol)) / ((-diff(5)+8*diff(4)-8*diff(2)+diff(1)) / (12*delta) / ucvol))
3691 
3692 end do
3693 
3694 end subroutine effective_potential_checkDEV

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

PARENTS

      m_effective_potential_file

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

202 subroutine effective_potential_init(crystal,eff_pot,energy,ifcs,ncoeff,nqpt,comm,&
203 &                                   coeffs,dynmat,elastic_constants,elastic3rd,&
204 &                                   elastic_displacement,epsilon_inf,&
205 &                                   fcart,strain_coupling,strten,name,phonon_strain,&
206 &                                   polynomial_conf,phfrq,qpoints,has_anharmonicsTerms,&
207 &                                   supercell,zeff)
208 
209 
210 !This section has been created automatically by the script Abilint (TD).
211 !Do not modify the following lines by hand.
212 #undef ABI_FUNC
213 #define ABI_FUNC 'effective_potential_init'
214 !End of the abilint section
215 
216  implicit none
217 
218 !Arguments ------------------------------------
219 !scalars
220  integer,intent(in) :: comm,ncoeff,nqpt
221  real(dp),intent(in):: energy
222  logical,optional,intent(in) :: has_anharmonicsTerms
223 !arrays
224  type(crystal_t),intent(in) :: crystal
225  type(effective_potential_type), intent(out) :: eff_pot
226  type(ifc_type),intent(in) :: ifcs
227  character(len=fnlen), optional,intent(in) :: name
228  real(dp),optional,intent(in) :: epsilon_inf(3,3),elastic_constants(6,6)
229  real(dp),optional,intent(in) :: dynmat(:,:,:,:,:,:),qpoints(:,:),phfrq(:,:)
230  real(dp),optional,intent(in) :: strain_coupling(:,:,:),zeff(:,:,:)
231  type(supercell_type),optional,intent(in) :: supercell
232  type(ifc_type),optional,intent(in) :: phonon_strain(6)
233  real(dp),optional,intent(in) :: elastic3rd(6,6,6)
234  real(dp),optional,intent(in) :: elastic_displacement(6,6,3,crystal%natom)
235  type(polynomial_coeff_type),optional :: coeffs(ncoeff)
236  real(dp),optional,intent(in) :: strten(6)
237  real(dp),optional,intent(in) :: fcart(3,crystal%natom)
238  type(polynomial_conf_type),optional,intent(in) :: polynomial_conf
239 
240 !Local variables-------------------------------
241 !scalar
242  character(len=500) :: msg
243 !arrays
244 
245 ! *************************************************************************
246 
247 !1-Free the effective potential before filling it
248  call effective_potential_free(eff_pot)
249 
250  if (present(name)) then
251    eff_pot%name = name
252  else
253    eff_pot%name = ''
254  end if
255 
256 !2-Perform some checks
257  if (crystal%natom < 1) then
258    write(msg, '(a,a,a,i10,a)' )&
259 &   'The cell must have at least one atom.',ch10,&
260 &   'The number of atom is  ',crystal%natom,'.'
261    MSG_BUG(msg)
262  end if
263 
264  if (crystal%ntypat < 1) then
265    write(msg, '(a,a,a,i10,a)' )&
266 &   'The cell must have at least one type of atom.',ch10,&
267 &   'The number of type of atom is  ',crystal%ntypat,'.'
268    MSG_BUG(msg)
269  end if
270 
271 !3-Fill energy of the crystal (hartree)
272  eff_pot%energy = energy
273 
274 !1-Fill the crystal
275 !Warning znucl is dimension with ntypat = nspsp hence alchemy is not supported here
276  call crystal_init(crystal%amu,eff_pot%crystal,crystal%space_group,crystal%natom,&
277 &                  crystal%npsp,crystal%ntypat,crystal%nsym,crystal%rprimd,&
278 &                  crystal%typat,crystal%xred,crystal%zion,crystal%znucl,&
279 &                  crystal%timrev,.FALSE.,.FALSE.,crystal%title,&
280 &                  symrel=crystal%symrel,tnons=crystal%tnons,symafm=crystal%symafm)
281 
282 !4-Fill harmonic part
283  call harmonics_terms_init(eff_pot%harmonics_terms,ifcs,crystal%natom,ifcs%nrpt)
284 
285 !5-Init the anharmonics_terms to set the flag to false
286  call anharmonics_terms_init(eff_pot%anharmonics_terms,crystal%natom,ncoeff)
287 
288 !5-Fill optional inputs
289  ABI_ALLOCATE(eff_pot%fcart,(3,eff_pot%crystal%natom))
290  eff_pot%fcart = zero
291  if(present(fcart))then
292    eff_pot%fcart = fcart
293  end if
294 
295  if(present(elastic_constants))then
296    eff_pot%harmonics_terms%elastic_constants(:,:) = elastic_constants
297  end if
298 
299  if(present(epsilon_inf))then
300    eff_pot%harmonics_terms%epsilon_inf(:,:) = epsilon_inf(:,:)
301  end if
302 
303  if(present(dynmat).and.present(qpoints).and.present(phfrq))then
304    call harmonics_terms_setDynmat(dynmat,eff_pot%harmonics_terms,crystal%natom,nqpt,phfrq,qpoints)
305  end if
306 
307  if(present(strain_coupling))then
308    call harmonics_terms_setInternalStrain(eff_pot%harmonics_terms,crystal%natom,strain_coupling)
309  end if
310 
311  if(present(zeff))then
312    call harmonics_terms_setEffectiveCharges(eff_pot%harmonics_terms,crystal%natom,zeff)
313  end if
314 
315  eff_pot%strten = zero
316  if(present(strten))then
317    eff_pot%strten(:) = strten(:)
318  end if
319 
320 !Set the flag for the strain coupling
321  if(present(has_anharmonicsTerms)) then
322    eff_pot%has_anharmonicsTerms = has_anharmonicsTerms
323  else
324    eff_pot%has_anharmonicsTerms = .false.
325  end if
326 
327 !Allocation of phonon strain coupling array (3rd order)
328  if(present(phonon_strain).and.has_anharmonicsTerms) then
329    call anharmonics_terms_setStrainPhononCoupling(eff_pot%anharmonics_terms,crystal%natom,&
330 &                                                 phonon_strain)
331  end if
332 
333 !Set the 3rd order elastic tensor
334  if(present(elastic3rd).and.has_anharmonicsTerms)then
335    call anharmonics_terms_setElastic3rd(eff_pot%anharmonics_terms,elastic3rd)
336  end if
337 
338 !Allocation of 3rd order with respecto to 2 strain and 1 atomic displacement
339  if(present(elastic_displacement).and.has_anharmonicsTerms)then
340    call anharmonics_terms_setElasticDispCoupling(eff_pot%anharmonics_terms,crystal%natom,&
341 &                                                elastic_displacement)
342  end if
343 
344 !Allocation of the coefficients
345  if(present(coeffs))then
346    if(ncoeff /= size(coeffs))then
347      write(msg, '(a)' )&
348 &        ' ncoeff has not the same size than coeffs array, '
349      MSG_BUG(msg)
350    end if
351    call effective_potential_setCoeffs(coeffs,eff_pot,ncoeff)
352  end if
353 
354  if(present(supercell))then
355    call effective_potential_setSupercell(eff_pot,comm,supercell=supercell)
356  else
357 !   call init_supercell(eff_pot%crystal%natom, (/1,0,0, 0,1,0,  0,0,1/), eff_pot%crystal%rprimd,&
358 !&                      eff_pot%crystal%typat, eff_pot%crystal%xcart, eff_pot%crystal%znucl, eff_pot%supercell)
359    call effective_potential_setSupercell(eff_pot,comm,ncell=(/1,1,1/))
360  end if
361 
362 !Set the confinement potential
363  if(present(polynomial_conf)) then
364    call effective_potential_setConfinement(polynomial_conf%cutoff_disp,polynomial_conf%cutoff_strain,&
365 &                                          eff_pot,polynomial_conf%factor_disp,&
366 &                                          polynomial_conf%factor_strain,polynomial_conf%ndisp,&
367 &                                          polynomial_conf%power_disp,polynomial_conf%power_strain,&
368 &                                          polynomial_conf%need_confinement)
369  end if
370 
371 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 informations 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

PARENTS

      m_effective_potential

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

403 subroutine effective_potential_initmpi(eff_pot,comm)
404 
405 
406 !This section has been created automatically by the script Abilint (TD).
407 !Do not modify the following lines by hand.
408 #undef ABI_FUNC
409 #define ABI_FUNC 'effective_potential_initmpi'
410 !End of the abilint section
411 
412  implicit none
413 
414 !Arguments ------------------------------------
415 !scalars
416  type(effective_potential_type),intent(inout)  :: eff_pot
417  integer,intent(in) :: comm
418 !arrays
419 
420 !Local variables-------------------------------
421 !scalars
422  integer :: ndiv
423  integer :: ncell
424 !array
425  integer :: cell_number(3)
426  character(len=500) :: msg
427 ! ***********************************************************************
428 
429 !Set the number of cell in the supercell
430  cell_number(1) = eff_pot%supercell%rlatt(1,1)
431  cell_number(2) = eff_pot%supercell%rlatt(2,2)
432  cell_number(3) = eff_pot%supercell%rlatt(3,3)
433  ncell = product(cell_number(:))
434 
435 !Do some checks
436  if (any(cell_number <= 0).or.ncell<=0) then
437    write(msg,'(a,a)')' No supercell found for setting'
438    MSG_ERROR(msg)
439  end if
440 
441 !First mpi_ifc
442  ndiv = 1
443  call effpot_mpi_free(eff_pot%mpi_ifc)
444  call effpot_mpi_init(eff_pot%harmonics_terms%ifcs%cell,cell_number,eff_pot%mpi_ifc,&
445 &                     eff_pot%crystal%natom,ndiv,eff_pot%harmonics_terms%ifcs%nrpt,comm)
446 
447 !Second mpi_coeff
448  ndiv = 1
449  call effpot_mpi_free(eff_pot%mpi_coeff)
450  call effpot_mpi_init(eff_pot%harmonics_terms%ifcs%cell,cell_number,eff_pot%mpi_coeff,&
451 &                     eff_pot%crystal%natom,ndiv,eff_pot%harmonics_terms%ifcs%nrpt,comm)
452 
453 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

PARENTS

      m_effective_potential,multibinit

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

1365 subroutine effective_potential_setConfinement(cutoff_disp,cutoff_strain,eff_pot,factor_disp,&
1366 &                                             factor_strain,ndisp,power_disp,power_strain,&
1367 &                                             need_confinement)
1368 
1369 
1370 !This section has been created automatically by the script Abilint (TD).
1371 !Do not modify the following lines by hand.
1372 #undef ABI_FUNC
1373 #define ABI_FUNC 'effective_potential_setConfinement'
1374 !End of the abilint section
1375 
1376  implicit none
1377 
1378 !Arguments ------------------------------------
1379 !scalars
1380  integer, intent(in) :: power_disp,power_strain,ndisp
1381  real(dp),intent(in) :: factor_disp,factor_strain
1382  logical,optional,intent(in)  :: need_confinement
1383 !arrays
1384  real(dp),intent(in) :: cutoff_disp(ndisp),cutoff_strain(6)
1385  type(effective_potential_type),intent(inout) :: eff_pot
1386 !Local variables-------------------------------
1387 !scalar
1388  logical  :: need_confinement_tmp = .FALSE.
1389 !arrays
1390  character(len=500) :: msg
1391 
1392 ! *************************************************************************
1393 
1394 !Checks
1395  if (ndisp <= 0) then
1396    write(msg,'(a,a)')' ndisp can not be inferior or equal to zero'
1397    MSG_ERROR(msg)
1398  end if
1399 
1400 !First free the type
1401  call  polynomial_conf_free(eff_pot%confinement)
1402 
1403  if (present(need_confinement)) need_confinement_tmp = need_confinement
1404 
1405  call  polynomial_conf_init(cutoff_disp,cutoff_strain,factor_disp,factor_strain,ndisp,&
1406 &                           eff_pot%confinement,power_disp,power_strain,&
1407 &                           need_confinement=need_confinement_tmp)
1408 
1409 
1410 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

PARENTS

      m_effective_potential,m_effective_potential_file,mover_effpot

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

1438 subroutine effective_potential_setSupercell(eff_pot,comm,ncell,supercell)
1439 
1440 
1441 !This section has been created automatically by the script Abilint (TD).
1442 !Do not modify the following lines by hand.
1443 #undef ABI_FUNC
1444 #define ABI_FUNC 'effective_potential_setSupercell'
1445 !End of the abilint section
1446 
1447  implicit none
1448 
1449 !Arguments ------------------------------------
1450 !scalars
1451  integer,intent(in) :: comm
1452 !arrays
1453  type(effective_potential_type),intent(inout) :: eff_pot
1454  integer,optional,intent(in) :: ncell(3)
1455  type(supercell_type),optional,intent(in) :: supercell
1456 !Local variables-------------------------------
1457 !scalar
1458 !arrays
1459  character(len=500) :: msg
1460 
1461 ! *************************************************************************
1462 
1463 !Checks
1464  if (.not.present(supercell).and..not.present(ncell)) then
1465    write(msg,'(a,a)')' You should at least set ncell of supercell type'
1466    MSG_ERROR(msg)
1467  end if
1468 
1469  call destroy_supercell(eff_pot%supercell)
1470 
1471  if(present(supercell))then
1472    call copy_supercell(supercell,eff_pot%supercell)
1473  else
1474    call init_supercell(eff_pot%crystal%natom, (/ncell(1),0,0,  0,ncell(2),0,  0,0,ncell(3)/), &
1475 &                      eff_pot%crystal%rprimd,eff_pot%crystal%typat,eff_pot%crystal%xcart,&
1476 &                      eff_pot%crystal%znucl, eff_pot%supercell)
1477  end if
1478 
1479 !Initialisation of new mpi over supercell
1480  call effective_potential_initmpi(eff_pot,comm)
1481 
1482 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

105  type, public :: effective_potential_type
106 
107    character(len=fnlen) :: name
108 !     Name of the molecule (CaTiO3,...)
109 
110    type(crystal_t) :: crystal
111 !    crystal type
112 !    contains all information of the crystal
113 
114    real(dp):: energy
115 !     Energy of the system (Hatree)
116 
117    real(dp), allocatable :: fcart(:,:)
118 !    forces(3,natom)
119 !    initial cartesian forces of the system
120 
121    real(dp) :: strten(6)
122 !    strten(6)
123 !    initial stresses (Ha/bohr^3) of the system
124 
125    type(harmonics_terms_type) :: harmonics_terms
126 !     type with all information for harmonics terms
127 
128    type(anharmonics_terms_type) :: anharmonics_terms
129 !     type with all information for anharmonics terms
130 
131    type(polynomial_conf_type) :: confinement
132 !     type with all the informations for the confinement
133 
134    type(supercell_type) :: supercell
135 !     super cell type
136 !     Store all the information of the suppercell
137 
138    logical :: has_anharmonicsTerms
139 !     True : the aharmonic part is present
140 
141 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
142 ! This is for the parallelisation over the supercell
143    type(effpot_mpi_type) :: mpi_ifc
144 !  effpot_mpi_type with all the information for the IFC paralellisation
145 
146    type(effpot_mpi_type) :: mpi_coeff
147 !  effpot_mpi_type with all the information for the polynomial coefficients paralellisation
148 
149  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

COPYRIGHT

 Copyright (C) 2000-2018 ABINIT group (AM)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

PARENTS

      compute_anharmonics

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

2144 subroutine effective_potential_writeAbiInput(eff_pot,filename,strain)
2145 
2146 
2147 !This section has been created automatically by the script Abilint (TD).
2148 !Do not modify the following lines by hand.
2149 #undef ABI_FUNC
2150 #define ABI_FUNC 'effective_potential_writeAbiInput'
2151 !End of the abilint section
2152 
2153   implicit none
2154 
2155 !Arguments ------------------------------------
2156 !scalars
2157   type(strain_type),optional,intent(in) :: strain
2158   character(len=fnlen),optional,intent(in) :: filename
2159 !arrays
2160   type(effective_potential_type), intent(in) :: eff_pot
2161 !Local variables-------------------------------
2162 !scalar
2163  integer :: unit = 20
2164  character(len=500) :: msg
2165  character(len=fnlen) :: namefile
2166 !arrays
2167  real(dp) :: xred(3,eff_pot%crystal%natom)
2168  type(strain_type) :: strain_tmp
2169 
2170 ! ************************************************************************
2171 
2172  if(present(strain)) then
2173    strain_tmp = strain
2174  else
2175    call strain_init(strain_tmp)
2176  end if
2177 
2178 ! try to open the file
2179  if(present(filename)) then
2180    namefile=filename
2181  else
2182    if(eff_pot%name /='') then
2183      write(namefile,'(a)') 'structure_'//trim(eff_pot%name)
2184    else
2185      write(namefile,'(a)') 'structure'
2186    end if
2187    if (strain_tmp%name/='') then
2188      write(namefile,'(a,a,a,a,a,a,a)') trim(namefile)//"_"//trim(strain_tmp%name)//"_"//&
2189 &        trim(itoa(strain_tmp%direction)),"_"//trim(ftoa(strain_tmp%delta))
2190    end if
2191 
2192    namefile=trim(namefile)//".in"
2193 
2194  end if
2195 
2196  call isfile(namefile,'new')
2197 
2198  if (open_file(namefile,msg,unit=unit,form="formatted",status="new",action="write") /= 0) then
2199    MSG_ERROR(msg)
2200  end if
2201 
2202   write(msg,'(a,a,a,a)')ch10,&
2203  &   ' Generation of the input file in ',namefile,ch10
2204   call wrtout(ab_out,msg,'COLL')
2205   call wrtout(std_out,msg,'COLL')
2206 
2207   write(unit,'("#Abinit Input for DFPT, this file contrains the keyword")')
2208   write(unit,'("#To run DFPT calculation of")')
2209   write(unit,'(a)') trim(eff_pot%name)
2210   if (strain_tmp%direction /= 0) then
2211     write(unit,'("# With a perturbation  ")',advance="no")
2212     write(unit,'(a)',advance="no") trim(strain_tmp%name)
2213     write(unit,'(" in the direction : ")',advance="no")
2214     write(unit,'(a)',advance='no') trim(itoa(strain_tmp%direction))
2215     write(unit,'(" with the deformation : ")',advance="no")
2216     write(unit,'(a)') trim(ftoa(strain_tmp%delta))
2217   end if
2218 
2219   write(unit,'("")')
2220   write(unit,'("ndtset 1 jdtset 1 2 3")')
2221   write(unit,'("")')
2222   write(unit,'("#DATASET1 GROUND  STATE")')
2223   write(unit,'("tolwfr1 = 1d-15")')
2224   write(unit,'(" prtwf1 = 1")')
2225   write(unit,'(" nline1 = 5")')
2226 
2227   write(unit,'("")')
2228   write(unit,'("#DATASET2 DDK PERTURBATION")')
2229   write(unit,'("getwfk2 =  1")')
2230   write(unit,'("  iscf2 = -3")')
2231   write(unit,'(" nline2 =  15")')
2232   write(unit,'("nnsclo2 =  5")')
2233   write(unit,'("kptopt2 =  2")')
2234   write(unit,'("  nqpt2 =  1")')
2235   write(unit,'("   qpt2 =  0 0 0 ")')
2236   write(unit,'("rfelfd2 =  2")')
2237   write(unit,'(" rfdir2 =  1 1 1 ")')
2238   write(unit,'("tolwfr2 =  1.0d-20 ")')
2239   write(unit,'(" prtwf2 =  1 ")')
2240   write(unit,'(" ")')
2241 
2242   write(unit,'("#DATASET3 RF")')
2243   write(unit,'(" getddk3 =  2")')
2244   write(unit,'(" getwfk3 =  1")')
2245   write(unit,'("   iscf3 =  7")')
2246   write(unit,'(" kptopt3 =  2")')
2247   write(unit,'("   nqpt3 =  1")')
2248   write(unit,'("    qpt3 =  0 0 0")')
2249   write(unit,'(" rfphon3 =  1")')
2250   write(unit,'("rfatpol3 =  1 ")',advance='no')
2251   write(unit,'(a)') itoa(eff_pot%crystal%natom)
2252   write(unit,'(" rfelfd3 =  3")')
2253   write(unit,'(" rfstrs3 =  3")')
2254   write(unit,'("  rfdir3 =  1 1 1")')
2255   write(unit,'(" tolvrs3 =  1.0d-8")')
2256   write(unit,'("")')
2257 
2258   write(unit,'("#STRUCTURE")')
2259   write(unit,'(" natom = ")',advance='no')
2260   write(unit,'(a)') itoa(eff_pot%crystal%natom)
2261   write(unit,'(" znucl =")',advance='no')
2262   write(unit,'(10(F4.0))') (eff_pot%crystal%znucl)
2263   write(unit,'("ntypat = ")',advance='no')
2264   write(unit,'(a)') itoa(eff_pot%crystal%ntypat)
2265   write(unit,'(" typat = ")',advance='no')
2266   write(unit,'(10(I2))') (eff_pot%crystal%typat)
2267   write(unit,'(" acell = 1 1 1")')
2268   write(unit,'(" rprim  ")')
2269   write(unit,'(3(F20.10))') (matmul(eff_pot%crystal%rprimd,strain%strain))
2270   write(unit,'("  xred  ")')
2271   call xcart2xred(eff_pot%crystal%natom,eff_pot%crystal%rprimd,eff_pot%crystal%xcart,xred)
2272   write(unit,'(3(F15.10))') (xred)
2273   write(unit,'(" ")')
2274 
2275   write(unit,'("#SCF")')
2276   write(unit,'("     ecut = ")')
2277   write(unit,'("pawecutdg = ")')
2278   write(unit,'("   ecutsm = ")')
2279   write(unit,'("   tolvrs = ")')
2280   write(unit,'("    nband = ")')
2281   write(unit,'("      ixc = ")')
2282   write(unit,'("   occopt = ")')
2283   write(unit,'("    nstep = ")')
2284   write(unit,'("   kptopt = ")')
2285   write(unit,'("    ngkpt =  ")')
2286   write(unit,'("")')
2287   write(unit,'("    prtwf 0 prtden 0 prtdos 0")')
2288 
2289   close(unit)
2290 
2291 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

COPYRIGHT

 Copyright (C) 2000-2018 ABINIT group (AM)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

PARENTS

      multibinit

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

3728 subroutine effective_potential_writeNETCDF(eff_pot,option,filename)
3729 
3730 
3731 !This section has been created automatically by the script Abilint (TD).
3732 !Do not modify the following lines by hand.
3733 #undef ABI_FUNC
3734 #define ABI_FUNC 'effective_potential_writeNETCDF'
3735 !End of the abilint section
3736 
3737   implicit none
3738 
3739 !Arguments ------------------------------------
3740 !scalars
3741   integer, intent(in) :: option
3742   character(len=fnlen),optional,intent(in) :: filename
3743 !arrays
3744   type(effective_potential_type), intent(in) :: eff_pot
3745 
3746 !Local variables-------------------------------
3747 !scalar
3748  integer :: amu_id,bec_id,ifccell_id,epsinf_id,elastic_id
3749  integer :: ifc_id,ifcs_id,natom_id,ntypat_id,nrpt_id,npsp_id,typat_id
3750  integer :: six_id,two_id,xyz_id,znucl_id
3751  integer :: ncerr,ncid,npsp
3752  integer :: dimCids(2),dimEids(2),dimIids(6),dimPids(1),dimRids(2),dimXids(2)
3753  integer :: etotal_id,rprimd_id,xcart_id
3754  character(len=500) :: msg
3755  character(len=fnlen) :: namefile
3756 !arrays
3757  real(dp) :: strain(9,6)
3758 
3759 ! *************************************************************************
3760 
3761 #if defined HAVE_NETCDF
3762 
3763  strain(:,1) = (/1,0,0,0,0,0,0,0,0/)
3764  strain(:,2) = (/0,0,0,0,1,0,0,0,0/)
3765  strain(:,3) = (/0,0,0,0,0,0,0,0,1/)
3766  strain(:,4) = half*(/0,0,0,0,0,1,0,1,0/)
3767  strain(:,5) = half*(/0,0,1,0,0,0,1,0,0/)
3768  strain(:,6) = half*(/0,1,0,1,0,0,0,0,0/)
3769 
3770 !Print only the reference system in xml format
3771  if (option == 1) then
3772 
3773    if(present(filename)) then
3774      namefile=filename
3775    else
3776      namefile='ref.nc'
3777    end if
3778 
3779    call isfile(namefile,'new')
3780 
3781    write(msg,'(a,a,a)')ch10,&
3782 &   ' Generation of the xml file for the reference structure in ',namefile
3783 
3784    call wrtout(ab_out,msg,'COLL')
3785    call wrtout(std_out,msg,'COLL')
3786 
3787 !  1. Create netCDF file
3788    ncerr = nf90_create(path=trim(namefile),cmode=NF90_CLOBBER, ncid=ncid)
3789    NCF_CHECK_MSG(ncerr,"create netcdf history file")
3790 
3791 !  2. Define dimensions
3792    ncerr = nf90_def_dim(ncid,"natom",eff_pot%crystal%natom,natom_id)
3793    NCF_CHECK_MSG(ncerr," define dimension natom")
3794 
3795    ncerr = nf90_def_dim(ncid,"ntypat",eff_pot%crystal%ntypat,ntypat_id)
3796    NCF_CHECK_MSG(ncerr," define dimension ntypat")
3797 
3798    ncerr = nf90_def_dim(ncid,"nrpt",eff_pot%harmonics_terms%ifcs%nrpt,nrpt_id)
3799    NCF_CHECK_MSG(ncerr," define dimension ntypat")
3800 
3801    ncerr = nf90_def_var(ncid, "typat", NF90_INT, natom_id, typat_id)
3802    NCF_CHECK_MSG(ncerr," define variable typat")
3803 
3804    npsp = size(eff_pot%crystal%znucl)
3805    if (npsp /= eff_pot%crystal%ntypat) then
3806      MSG_WARNING("HIST file does not support alchemical mixing!")
3807    end if
3808    ncerr = nf90_def_dim(ncid,"npsp",npsp,npsp_id)
3809    NCF_CHECK_MSG(ncerr," define dimension npsp")
3810 
3811    ncerr = nf90_def_var(ncid, "znucl", NF90_DOUBLE, npsp_id, znucl_id)
3812    NCF_CHECK_MSG(ncerr," define variable znucl")
3813 
3814    ncerr = nf90_def_dim(ncid,"xyz",3,xyz_id)
3815    NCF_CHECK_MSG(ncerr," define dimension xyz")
3816 
3817    ncerr = nf90_def_dim(ncid,"six",6,six_id)
3818    NCF_CHECK_MSG(ncerr," define dimension six")
3819 
3820    ncerr = nf90_def_dim(ncid,"two",2,two_id)
3821    NCF_CHECK_MSG(ncerr," define dimension two")
3822 
3823 !  Dimensions for xcart,xred,fcart,fred and vel
3824    dimXids = (/ xyz_id, natom_id /)
3825 !  Dimensions for rprimd
3826    dimRids = (/ xyz_id, xyz_id /)
3827 !  Dimensions for ifc
3828    dimIids = (/ 2, xyz_id, natom_id, xyz_id, natom_id, nrpt_id /)
3829 !  Dimensions for position
3830    dimPids = (/ xyz_id /)
3831 !  Dimension for elastic constant
3832    dimEids = (/six_id,six_id/)
3833 !  Dimension for cell
3834    dimCids = (/nrpt_id,3/)
3835 
3836 !  3. Define variables and their attributes (units and mnemonics)
3837    call ab_define_var(ncid, (/1/), etotal_id, NF90_DOUBLE,&
3838 &   "energy","Energy of the reference structure","Ha" )
3839 
3840    call ab_define_var(ncid, dimRids, rprimd_id, NF90_DOUBLE,&
3841 &   "rprimd","Real space PRIMitive translations, Dimensional","bohr" )
3842 
3843    call ab_define_var(ncid, dimRids, epsinf_id, NF90_DOUBLE,&
3844 &   "epsilon_inf","Dielectric tensor, Dimensional","epsilon_inf" )
3845 
3846    call ab_define_var(ncid, dimEids, elastic_id, NF90_DOUBLE,&
3847 &   "elastic","Elastic Constants, Dimensional","Ha" )
3848 
3849    call ab_define_var(ncid, dimRids, bec_id, NF90_DOUBLE,&
3850 &   "bec","Born Effective Charges, Dimensional","abs(e)" )
3851 
3852    call ab_define_var(ncid, dimXids, xcart_id, NF90_DOUBLE,&
3853 &   "xcart","vectors (X) of atom positions in CARTesian coordinates","bohr" )
3854 
3855    call ab_define_var(ncid, dimIids, ifcs_id, NF90_DOUBLE,&
3856 &   "IFCs","Interatomic Forces Constantes in real spaces (short range), Dimensional","Hatree/bohr**2" )
3857 
3858    call ab_define_var(ncid, dimIids, ifc_id, NF90_DOUBLE,&
3859 &   "IFC","Interatomic Forces Constantes in real spaces (total range), Dimensional","Hatree/bohr**2" )
3860 
3861    call ab_define_var(ncid, dimCids, ifccell_id, NF90_DOUBLE,&
3862 &   "cell","cell for the ifc, Dimensional","Dimensionless" )
3863 
3864    call ab_define_var(ncid, [ntypat_id], amu_id, NF90_DOUBLE,&
3865 &   "amu","Masses of each type of atom in atomic mass units", "" )
3866 
3867 !  4. End define mode
3868    ncerr = nf90_enddef(ncid)
3869    NCF_CHECK_MSG(ncerr," end define mode")
3870 
3871 !  5. Write variables
3872    ncerr = nf90_put_var(ncid,etotal_id, eff_pot%energy)
3873    NCF_CHECK_MSG(ncerr," write variable energy")
3874 
3875    ncerr = nf90_put_var(ncid,rprimd_id, eff_pot%crystal%rprimd)
3876    NCF_CHECK_MSG(ncerr," write variable rprimd")
3877 
3878    ncerr = nf90_put_var(ncid,epsinf_id, eff_pot%harmonics_terms%epsilon_inf)
3879    NCF_CHECK_MSG(ncerr," write variable epsilon_inf")
3880 
3881    ncerr = nf90_put_var(ncid,elastic_id , eff_pot%harmonics_terms%elastic_constants)
3882    NCF_CHECK_MSG(ncerr," write variable elastic_constant")
3883 
3884    ncerr = nf90_put_var(ncid, bec_id, eff_pot%harmonics_terms%zeff)
3885    NCF_CHECK_MSG(ncerr," write variable bec")
3886 
3887    ncerr = nf90_put_var(ncid,xcart_id, eff_pot%crystal%xcart)
3888    NCF_CHECK_MSG(ncerr," write variable xcart")
3889 
3890    ncerr = nf90_put_var(ncid,ifccell_id, eff_pot%harmonics_terms%ifcs%cell)
3891    NCF_CHECK_MSG(ncerr," write variable cell")
3892 
3893    ncerr = nf90_put_var(ncid,ifcs_id, eff_pot%harmonics_terms%ifcs%short_atmfrc)
3894    NCF_CHECK_MSG(ncerr," write variable short ifc")
3895 
3896    ncerr = nf90_put_var(ncid,ifc_id, eff_pot%harmonics_terms%ifcs%atmfrc)
3897    NCF_CHECK_MSG(ncerr," write variable total ifc")
3898 
3899 
3900 !  6. Close NetCDF file
3901    ncerr = nf90_close(ncid)
3902    NCF_CHECK_MSG(ncerr," close netcdf history file")
3903  end if
3904 
3905 #endif
3906 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

COPYRIGHT

 Copyright (C) 2000-2018 ABINIT group (AM)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

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

PARENTS

      multibinit

CHILDREN

      ab_define_var,isfile,wrtout

SOURCE

1804 subroutine effective_potential_writeXML(eff_pot,option,filename,prt_dipdip)
1805 
1806 
1807 !This section has been created automatically by the script Abilint (TD).
1808 !Do not modify the following lines by hand.
1809 #undef ABI_FUNC
1810 #define ABI_FUNC 'effective_potential_writeXML'
1811 !End of the abilint section
1812 
1813   implicit none
1814 
1815 !Arguments ------------------------------------
1816 !scalars
1817   integer, intent(in) :: option
1818   character(len=fnlen),optional,intent(in) :: filename
1819   logical,optional,intent(in) :: prt_dipdip
1820 !arrays
1821   type(effective_potential_type), intent(in) :: eff_pot
1822 
1823 !Local variables-------------------------------
1824 !scalar
1825  integer :: ii,ia,ib,jj
1826  integer :: iqpt,irpt,mu,nu
1827  integer :: unit_xml
1828  character(len=500) :: msg
1829  character(len=fnlen) :: namefile
1830  character(len=10) :: natom
1831  logical :: new_file = .FALSE.,need_prtdipdip = .FALSE.
1832 !arrays
1833  real(dp) :: strain(9,6)
1834 
1835 ! *************************************************************************
1836 
1837  strain(:,1) = (/1,0,0,0,0,0,0,0,0/)
1838  strain(:,2) = (/0,0,0,0,1,0,0,0,0/)
1839  strain(:,3) = (/0,0,0,0,0,0,0,0,1/)
1840  strain(:,4) = half*(/0,0,0,0,0,1,0,1,0/)
1841  strain(:,5) = half*(/0,0,1,0,0,0,1,0,0/)
1842  strain(:,6) = half*(/0,1,0,1,0,0,0,0,0/)
1843 
1844  unit_xml = get_unit()
1845  need_prtdipdip = .TRUE.
1846  if(present(prt_dipdip)) need_prtdipdip = prt_dipdip
1847 !Print only the reference system in xml format
1848  if (option ==  1 .or. option == 2 .or. option ==3) then
1849 
1850 !  convert natom in character
1851    write (natom,'(I9)') eff_pot%crystal%natom
1852 
1853 !  Compute the name of the XML file
1854    if(present(filename)) then
1855      namefile=replace(trim(filename),".out","")
1856      select case(option)
1857      case(1)
1858        namefile=trim(namefile)//"_model.xml"
1859      case(2)
1860        namefile=trim(namefile)//"_sys.xml"
1861      case(3)
1862        namefile=trim(namefile)//"_sys.xml"
1863      end select
1864    else
1865      namefile='system.xml'
1866    end if
1867 
1868    call isfile(namefile,'new')
1869 
1870    if (open_file(namefile,msg,unit=unit_xml,form="formatted",&
1871 &      status="new",action="write") /= 0) then
1872      MSG_ERROR(msg)
1873    end if
1874 
1875    write(msg,'(a,a,a)')ch10,&
1876  &       ' Generation of the xml file for the model in ',trim(namefile)
1877 
1878    call wrtout(ab_out,msg,'COLL')
1879    call wrtout(std_out,msg,'COLL')
1880 
1881 !  Write header
1882    WRITE(unit_xml,'("<?xml version=""1.0"" ?>")')
1883    WRITE(unit_xml,'("<System_definition>")')
1884 
1885    WRITE(unit_xml,'("  <energy>")')
1886    WRITE(unit_xml,'(E23.14)') (eff_pot%energy)
1887    WRITE(unit_xml,'("  </energy>")')
1888 
1889    WRITE(unit_xml,'("  <unit_cell units=""bohrradius"">")')
1890    do mu=1,3
1891      WRITE(unit_xml,'(3(E23.14))') (eff_pot%crystal%rprimd(mu,nu),nu=1,3)
1892    end do
1893    WRITE(unit_xml,'("  </unit_cell>")')
1894 
1895    WRITE(unit_xml,'("  <epsilon_inf units=""epsilon0"">")')
1896    WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%epsilon_inf)
1897    WRITE(unit_xml,'("  </epsilon_inf>")')
1898 
1899    WRITE(unit_xml,'("  <elastic units=""hartree"">")')
1900    WRITE(unit_xml,'(6(E23.14))') (eff_pot%harmonics_terms%elastic_constants)
1901    WRITE(unit_xml,'("  </elastic>")')
1902 
1903    do ia=1,eff_pot%crystal%natom
1904      WRITE(unit_xml,'("  <atom mass=""",1F10.5,""" massunits=""atomicmassunit"">")') &
1905 &      eff_pot%crystal%amu(eff_pot%crystal%typat(ia))
1906      WRITE(unit_xml,'("    <position units=""bohrradius"">")')
1907      WRITE(unit_xml,'(3(E23.14))') (eff_pot%crystal%xcart(:,ia))
1908      WRITE(unit_xml,'("    </position>")')
1909      WRITE(unit_xml,'("    <borncharge units=""abs(e)"">")')
1910      WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%zeff(:,:,ia))
1911      WRITE(unit_xml,'("    </borncharge>")')
1912      WRITE(unit_xml,'("  </atom>")')
1913    end do
1914 !  Print the Ifc short range for each cell data is array 3*natom*3*natom
1915 !  [ [x1 x2 ....]
1916 !    [y1 y2 ....] for atom 1
1917 !    [z1 z2 ....]
1918 !    [x1 x2 ....]
1919 !    [y1 y2 ....] for atom 2
1920 !    [z1 z2 ....]
1921 !    ....       ]
1922 ! Warning : The IFC are print in other order 1,mu,ia,nu,ib which is not fortran way
1923 !           We do like that to because the previous script was in python
1924 !           When you read ifc from XML file with fotran, you have to tranpose the matrix
1925 !
1926    do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt
1927      if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol20)) then
1928        WRITE(unit_xml,'("  <local_force_constant units=""hartree/bohrradius**2"">")')
1929        WRITE(unit_xml,'("    <data>")')
1930        do ia=1,eff_pot%crystal%natom
1931          do mu=1,3
1932            do ib=1,eff_pot%crystal%natom
1933              do  nu=1,3
1934                WRITE(unit_xml,'(e22.14)', advance="no")&
1935 &                         (eff_pot%harmonics_terms%ifcs%short_atmfrc(mu,ia,nu,ib,irpt))
1936              end do
1937            end do
1938            WRITE(unit_xml,'(a)')''
1939          end do
1940        end do
1941        WRITE(unit_xml,'("    </data>")')
1942        WRITE(unit_xml,'("    <cell>")')
1943        WRITE(unit_xml,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt))
1944        WRITE(unit_xml,'("    </cell>")')
1945        WRITE(unit_xml,'("  </local_force_constant>")')
1946      end if
1947 !    Print the IFC total for each cell, data is array 3*natom*3*natom
1948 !    [ [x1 x2 ....]
1949 !      [y1 y2 ....] for atom 1
1950 !      [z1 z2 ....]
1951 !      [x1 x2 ....]
1952 !      [y1 y2 ....] for atom 2
1953 !      [z1 z2 ....]
1954 !      ....       ]
1955      if(need_prtdipdip)then
1956        if(all(abs(eff_pot%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt))<tol20)) then
1957          if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol20)) then
1958            write(msg, '(a,a,a,a)' )&
1959 &         ' There is no total range but short range in your effective potential',ch10,&
1960 &         'Action: contact abinit group'
1961            MSG_BUG(msg)
1962          end if
1963        else
1964          WRITE(unit_xml,'("  <total_force_constant units=""hartree/bohrradius**2"">")')
1965          WRITE(unit_xml,'("    <data>")')
1966          do ia=1,eff_pot%crystal%natom
1967            do mu=1,3
1968              do ib=1,eff_pot%crystal%natom
1969                do nu=1,3
1970                  WRITE(unit_xml,'(e22.14)', advance="no")&
1971 &                   (eff_pot%harmonics_terms%ifcs%atmfrc(mu,ia,nu,ib,irpt))
1972                end do
1973              end do
1974              WRITE(unit_xml,'(a)')''
1975            end do
1976          end do
1977          WRITE(unit_xml,'("    </data>")')
1978          WRITE(unit_xml,'("    <cell>")')
1979          WRITE(unit_xml,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt))
1980          WRITE(unit_xml,'("    </cell>")')
1981          WRITE(unit_xml,'("    </total_force_constant>")')
1982        end if
1983      end if
1984    end do
1985 
1986    do iqpt=1,eff_pot%harmonics_terms%nqpt
1987      WRITE(unit_xml,'("  <phonon>")')
1988      WRITE(unit_xml,'("    <qpoint units=""2pi*G0"">")')
1989      WRITE(unit_xml,'(3(E23.14))') (eff_pot%harmonics_terms%qpoints(:,iqpt))
1990      WRITE(unit_xml,'("    </qpoint>")')
1991      WRITE(unit_xml,'("    <frequencies units=""reciprocal cm"">")')
1992      WRITE(unit_xml,'(3(e22.14))') (eff_pot%harmonics_terms%phfrq(:,iqpt))
1993      WRITE(unit_xml,'("    </frequencies>")')
1994      WRITE(unit_xml,'("    <dynamical_matrix units=""hartree/bohrradius**2"">")')
1995      do ia=1,eff_pot%crystal%natom
1996        do mu=1,3
1997          do ib = 1,eff_pot%crystal%natom
1998            do nu=1,3
1999              WRITE(unit_xml,'(e22.14)',advance='no')(eff_pot%harmonics_terms%dynmat(1,nu,ib,mu,ia,iqpt))
2000            end do
2001          end do
2002          WRITE(unit_xml,'(a)')''
2003        end do
2004      end do
2005      WRITE(unit_xml,'("    </dynamical_matrix>")')
2006      WRITE(unit_xml,'("  </phonon>")')
2007    end do
2008 ! if phonon/forces strain is computed
2009    jj = 1
2010    do ii = 1,6
2011      WRITE(unit_xml,'("  <strain_coupling voigt=""",I2,""">")') ii-1
2012      WRITE(unit_xml,'("    <strain>")')
2013      WRITE(unit_xml,'(6(e12.4))') (strain(:,jj))
2014      WRITE(unit_xml,'("    </strain>")')
2015      WRITE(unit_xml,'("    <correction_force units=""hartree/bohrradius"">")')
2016      do ia=1,eff_pot%crystal%natom
2017        do mu=1,3
2018          WRITE(unit_xml,'(e22.14)', advance="no")&
2019 &             (eff_pot%harmonics_terms%strain_coupling(ii,mu,ia))
2020        end do
2021        WRITE(unit_xml,'(a)')''
2022      end do
2023      WRITE(unit_xml,'("    </correction_force>")')
2024      if (eff_pot%has_anharmonicsTerms)then
2025        if (eff_pot%anharmonics_terms%has_elastic3rd) then
2026          WRITE(unit_xml,'("  <elastic3rd units=""hartree"">")')
2027          WRITE(unit_xml,'(6(E23.14))') (eff_pot%anharmonics_terms%elastic3rd(ii,:,:))
2028          WRITE(unit_xml,'("  </elastic3rd>")')
2029        end if
2030        if (eff_pot%anharmonics_terms%has_elastic_displ) then
2031          WRITE(unit_xml,'("    <correction_strain_force units=""hartree/bohrradius"">")')
2032          do ia=1,eff_pot%crystal%natom
2033            do mu=1,3
2034              do nu=1,6
2035                WRITE(unit_xml,'(e22.14)', advance="no")&
2036 &                   (eff_pot%anharmonics_terms%elastic_displacement(ii,nu,mu,ia))
2037              end do
2038            end do
2039            WRITE(unit_xml,'(a)')''
2040          end do
2041        end if
2042        WRITE(unit_xml,'("    </correction_strain_force>")')
2043        if (eff_pot%anharmonics_terms%has_strain_coupling) then
2044          do irpt=1,eff_pot%anharmonics_terms%phonon_strain(ii)%nrpt
2045            WRITE(unit_xml,'("    <correction_force_constant units=""hartree/bohrradius**2"">")')
2046            WRITE(unit_xml,'("      <data>")')
2047            do ia=1,eff_pot%crystal%natom
2048              do mu=1,3
2049                do ib=1,eff_pot%crystal%natom
2050                  do  nu=1,3
2051                    WRITE(unit_xml,'(e22.14)', advance="no")&
2052 &                      (eff_pot%anharmonics_terms%phonon_strain(ii)%atmfrc(mu,ia,nu,ib,irpt))
2053                  end do
2054                end do
2055                WRITE(unit_xml,'(a)')''
2056              end do
2057            end do
2058            WRITE(unit_xml,'("      </data>")')
2059            WRITE(unit_xml,'("      <cell>")')
2060            WRITE(unit_xml,'(3(I4))') (eff_pot%anharmonics_terms%phonon_strain(ii)%cell(:,irpt))
2061            WRITE(unit_xml,'("      </cell>")')
2062            WRITE(unit_xml,'("    </correction_force_constant>")')
2063          end do
2064        end if!End if has_straincouplitn
2065      end if!end Hasstrain_coupling
2066      WRITE(unit_xml,'("    </strain_coupling>")')
2067      jj = jj + 1
2068    end do!end mu
2069 
2070    if(option /=1)  WRITE(unit_xml,'("</System_definition>")')
2071 !  Close file
2072    CLOSE(unit_xml)
2073 
2074  end if!end option
2075 
2076 !Print the coefficients into XML file
2077  if (option==1 .or. option == 2 .or. option==4) then
2078 !   Compute the name of the XML file
2079    if(present(filename)) then
2080      namefile=replace(trim(filename),".out","")
2081      select case(option)
2082      case(1)
2083        new_file = .FALSE.
2084        namefile=trim(namefile)//"_model.xml"
2085      case(2)
2086        new_file = .TRUE.
2087        namefile=trim(namefile)//"_coeffs.xml"
2088      case(4)
2089        new_file = .TRUE.
2090        namefile=trim(namefile)//"_coeffs.xml"
2091      end select
2092    else
2093      namefile='coeffs.xml'
2094    end if
2095 
2096    if(eff_pot%anharmonics_terms%ncoeff > 0) then
2097      call polynomial_coeff_writeXML(eff_pot%anharmonics_terms%coefficients,&
2098 &                                   eff_pot%anharmonics_terms%ncoeff,namefile,unit=unit_xml,&
2099 &                                   newfile=new_file)
2100    end if
2101  end if!end option
2102 
2103  if(option==1)then
2104 !  add the end of the file in the case option==1
2105    open(unit=unit_xml,file=namefile,position="append")
2106    WRITE(unit_xml,'("</System_definition>")')
2107 !  Close file
2108    CLOSE(unit_xml)
2109  end if
2110 
2111 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

3065 pure function effective_potential_compare(e1,e2) result (res)
3066 !Arguments ------------------------------------
3067 
3068 !This section has been created automatically by the script Abilint (TD).
3069 !Do not modify the following lines by hand.
3070 #undef ABI_FUNC
3071 #define ABI_FUNC 'effective_potential_compare'
3072 !End of the abilint section
3073 
3074  implicit none
3075 
3076 !Arguments ------------------------------------
3077   type(effective_potential_type), intent(in) :: e1,e2
3078   logical :: res
3079 ! *************************************************************************
3080   res = .false.
3081   if(e1%crystal%natom==e2%crystal%natom.and.&
3082 &     e1%harmonics_terms%ifcs%nrpt==e2%harmonics_terms%ifcs%nrpt.and.&
3083 &     e1%crystal%ntypat==e2%crystal%ntypat.and.&
3084 &     e1%harmonics_terms%nqpt==e2%harmonics_terms%nqpt.and.&
3085 &     abs(e1%energy-e2%energy)<tol16.and.&
3086 &     abs(e1%crystal%ucvol-e2%crystal%ucvol)<tol16) then
3087     res = .true.
3088   end if
3089 
3090 end function effective_potential_compare