TABLE OF CONTENTS
- ABINIT/m_ddb
- m_ddb/asr_t
- m_ddb/asrq0_apply
- m_ddb/asrq0_free
- m_ddb/carteig2d
- m_ddb/carttransf
- m_ddb/chkin9
- m_ddb/ddb_bcast
- m_ddb/ddb_can_merge_blocks
- m_ddb/ddb_copy
- m_ddb/ddb_diagoq
- m_ddb/ddb_free
- m_ddb/ddb_from_file
- m_ddb/ddb_get_asrq0
- m_ddb/ddb_get_block
- m_ddb/ddb_get_d1matr
- m_ddb/ddb_get_d2eig
- m_ddb/ddb_get_d2matr
- m_ddb/ddb_get_d3matr
- m_ddb/ddb_get_dchidet
- m_ddb/ddb_get_dielt
- m_ddb/ddb_get_dielt_zeff
- m_ddb/ddb_get_etotal
- m_ddb/ddb_get_gred
- m_ddb/ddb_get_pel
- m_ddb/ddb_get_quadrupoles
- m_ddb/ddb_get_strten
- m_ddb/ddb_init
- m_ddb/ddb_lw_copy
- m_ddb/ddb_malloc
- m_ddb/ddb_malloc_d2eig
- m_ddb/ddb_merge_blocks
- m_ddb/ddb_read_block_txt
- m_ddb/ddb_read_d0E_nc
- m_ddb/ddb_read_d1E_nc
- m_ddb/ddb_read_d2E_nc
- m_ddb/ddb_read_d2eig
- m_ddb/ddb_read_d2eig_nc
- m_ddb/ddb_read_d2eig_txt
- m_ddb/ddb_read_d3E_nc
- m_ddb/ddb_read_nc
- m_ddb/ddb_read_txt
- m_ddb/ddb_set_brav
- m_ddb/ddb_set_d1matr
- m_ddb/ddb_set_d2eig
- m_ddb/ddb_set_d2eig_reshape
- m_ddb/ddb_set_d2matr
- m_ddb/ddb_set_d3matr
- m_ddb/ddb_set_etotal
- m_ddb/ddb_set_gred
- m_ddb/ddb_set_pel
- m_ddb/ddb_set_qpt
- m_ddb/ddb_set_strten
- m_ddb/ddb_set_typ
- m_ddb/ddb_symmetrize_and_transform
- m_ddb/ddb_to_dtset
- m_ddb/ddb_type
- m_ddb/ddb_write
- m_ddb/ddb_write_block_txt
- m_ddb/ddb_write_d2eig
- m_ddb/ddb_write_d2eig_nc
- m_ddb/ddb_write_d2eig_txt
- m_ddb/ddb_write_nc
- m_ddb/ddb_write_txt
- m_ddb/dtchi
- m_ddb/dtech9
- m_ddb/dtqdrp
- m_ddb/gamma9
- m_ddb/lwcart
- m_ddb/merge_ddb
- m_ddb/nlopt
- m_ddb/rdddb9
- m_ddb/symdm9
ABINIT/m_ddb [ Modules ]
NAME
m_ddb
FUNCTION
This module contains the declaration of data types and methods used to handle the blocks of data in DDB files: blkval, nrm, qpt, flg, and associated dimensions Main entry point for client code that needs to read the DDB data.
COPYRIGHT
Copyright (C) 2011-2024 ABINIT group (MJV, XG, MT, MM, MVeithen, MG, PB, JCC, SP, GA, MMignolet) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_ddb 26 27 use defs_basis 28 use m_abicore 29 use m_errors 30 use m_xmpi 31 use m_ddb_hdr 32 use m_dtset 33 use m_nctk 34 #ifdef HAVE_NETCDF 35 use netcdf 36 #endif 37 38 use m_io_tools, only : iomode_from_fname 39 use defs_datatypes, only : pseudopotential_type 40 use m_fstrings, only : sjoin, itoa, ktoa, endswith 41 use m_numeric_tools, only : mkherm 42 use m_symtk, only : mati3inv, matr3inv, littlegroup_q, symatm 43 use m_io_tools, only : get_unit 44 use m_copy, only : alloc_copy 45 use m_geometry, only : phdispl_cart2red, mkrdim, xred2xcart, metric 46 use m_crystal, only : crystal_t, crystal_init 47 use m_dynmat, only : cart29, d2sym3, cart39, d3sym, chneu9, asria_calc, asria_corr, asrprs, dfpt_phfrq, sytens 48 use m_pawtab, only : pawtab_type, pawtab_nullify, pawtab_free 49 use m_psps, only : psps_copy, psps_free 50 51 implicit none 52 53 private 54 55 public :: rdddb9 ! This routine reads the derivative database entirely, 56 public :: nlopt ! Output of all quantities related to third-order derivatives of the energy. 57 public :: chkin9 58 public :: carttransf ! Transform a second-derivative matrix (EIG2D) from reduced 59 ! coordinates to cartesian coordinates. 60 public :: lwcart ! Transform a 3rd order derivative tensor (long-wave) from reduced (actually 61 ! mixed since strain derivatives are already in cartesian) to cartesian 62 ! coordinates 63 public :: ddb_lw_copy ! Copy the ddb object after reading the long wave 3rd order derivatives 64 ! into a new ddb_lw and resizes ddb as for 2nd order derivatives 65 66 public :: symdm9 67 68 real(dp),public,parameter :: DDB_QTOL=2.0d-8 69 ! Tolerance for the identification of two wavevectors
m_ddb/asr_t [ Types ]
NAME
asr_t
FUNCTION
Object used to enforce the acoustic sum rule from the Dynamical matrix at Gamma. Wraps several approaches that can be activated via the `asr` option.
SOURCE
365 type,public :: asrq0_t 366 367 integer :: iblok = 0 368 ! Index of the Gamma block in the DDB. 369 ! Set to 0 if no block was found. Client code can use this flag to understand 370 ! if ASR can be enforced. 371 372 integer :: asr 373 ! Option for the application of the ASR (input variable). 374 375 integer :: natom 376 ! Number of atoms. 377 378 real(dp),allocatable :: d2asr(:,:,:,:,:) 379 ! d2asr,(2,3,natom,3,natom)) 380 ! In case the interatomic forces are not calculated, the 381 ! ASR-correction (d2asr) has to be determined here from the Dynamical matrix at Gamma. 382 383 ! singular, uinvers and vtinvers are allocated and used only if asr in [3,4] 384 ! i.e. Rotational invariance for 1D and 0D systems. dims=3*natom*(3*natom-1)/2 385 real(dp),allocatable :: singular(:) 386 ! singular,(1:dims)) 387 388 real(dp),allocatable :: uinvers(:,:) 389 ! uinvers,(1:dims,1:dims)) 390 391 real(dp),allocatable :: vtinvers(:,:) 392 ! vtinvers,(1:dims,1:dims)) 393 394 contains 395 396 procedure :: apply => asrq0_apply 397 ! Impose the acoustic sum rule based on the q=0 block found in the DDB file. 398 399 procedure :: free => asrq0_free 400 ! Free memory 401 402 end type asrq0_t
m_ddb/asrq0_apply [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
asrq0_apply
FUNCTION
Impose the acoustic sum rule based on the q=0 block found in the DDB file.
INPUTS
asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file. natom=Number of atoms per unit cell. mpert=Maximum number of perturbation (reported in ddb%mpert) msize=Maximum size of array ddb%val xcart(3,natom)=Atomic positions in Cartesian coordinates
SIDE EFFECTS
d2cart=matrix of second derivatives of total energy, in cartesian coordinates Input: Values stored in ddb% Output: Changed to enforce ASR.
SOURCE
4481 subroutine asrq0_apply(asrq0, natom, mpert, msize, xcart, d2cart) 4482 4483 !Arguments ------------------------------- 4484 !scalars 4485 integer,intent(in) :: natom, msize, mpert 4486 class(asrq0_t),intent(inout) :: asrq0 4487 !arrays 4488 real(dp),intent(in) :: xcart(3,natom) 4489 real(dp),intent(inout) :: d2cart(2,msize) 4490 4491 ! ************************************************************************ 4492 4493 if (asrq0%asr /= 0 .and. asrq0%iblok == 0) then 4494 ABI_WARNING("asr != 0 but DDB file does not contain q=Gamma. D(q) cannot be corrected") 4495 return 4496 end if 4497 4498 select case (asrq0%asr) 4499 case (0) 4500 return 4501 case (1,2,5) 4502 call asria_corr(asrq0%asr, asrq0%d2asr, d2cart, mpert, natom) 4503 case (3,4) 4504 ! Impose acoustic sum rule plus rotational symmetry for 0D and 1D systems 4505 call asrprs(asrq0%asr,2,3,asrq0%uinvers,asrq0%vtinvers,asrq0%singular,d2cart,mpert,natom,xcart) 4506 case default 4507 ABI_ERROR(sjoin("Wrong value for asr:", itoa(asrq0%asr))) 4508 end select 4509 4510 end subroutine asrq0_apply
m_ddb/asrq0_free [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
asrq0_free
FUNCTION
Free dynamic memory
SOURCE
4524 subroutine asrq0_free(asrq0) 4525 4526 !Arguments ------------------------------- 4527 class(asrq0_t),intent(inout) :: asrq0 4528 4529 ! ************************************************************************ 4530 4531 ! real 4532 ABI_SFREE(asrq0%d2asr) 4533 ABI_SFREE(asrq0%singular) 4534 ABI_SFREE(asrq0%uinvers) 4535 ABI_SFREE(asrq0%vtinvers) 4536 4537 end subroutine asrq0_free
m_ddb/carteig2d [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
carteig2d
FUNCTION
Transform a second-derivative matrix (EIG2D) from reduced coordinates to cartesian coordinates
INPUTS
blkflg(3,mpert,3,mpert,nblok)=
( 1 if the element of the dynamical matrix has been calculated ;
0 otherwise )
blkval(2,3,mpert,3,mpert,nblok)=DDB values
gprimd(3,3)=basis vector in the reciprocal space
iblok=number of the blok that will be transformed
mpert =maximum number of ipert
natom=number of atom
nblok=number of blocks (dimension of blkflg and blkval)
rprimd(3,3)=basis vector in the real space
OUTPUT
carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
2DTE matrix has been calculated correctly ; 0 otherwise )
d2cart(2,3,mpert,3,mpert)=
dynamical matrix, effective charges, dielectric tensor,....
all in cartesian coordinates
SOURCE
3153 subroutine carteig2d(blkflg,blkval,carflg,d2cart,gprimd,iblok,mpert,natom,nblok,rprimd) 3154 3155 !Arguments ------------------------------- 3156 !scalars 3157 integer,intent(in) :: iblok,mpert,natom,nblok 3158 !arrays 3159 integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok) 3160 integer,intent(out) :: carflg(3,mpert,3,mpert) 3161 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),gprimd(3,3),rprimd(3,3) 3162 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert) 3163 3164 !Local variables ------------------------- 3165 !scalars 3166 integer :: idir1,idir2,ii,ipert1,ipert2 3167 !arrays 3168 integer :: flg1(3),flg2(3) 3169 real(dp) :: vec1(3),vec2(3) 3170 3171 ! ********************************************************************* 3172 3173 ! First, copy the data blok in place. 3174 d2cart(:,:,:,:,:)=blkval(:,:,:,:,:,iblok) 3175 3176 ! Cartesian coordinates transformation (in two steps) 3177 ! First step 3178 do ipert1=1,mpert 3179 do ipert2=1,mpert 3180 do ii=1,2 3181 do idir1=1,3 3182 do idir2=1,3 3183 vec1(idir2)=d2cart(ii,idir1,ipert1,idir2,ipert2) 3184 ! Note here blkflg 3185 flg1(idir2)=blkflg(idir1,ipert1,idir2,ipert2,iblok) 3186 end do 3187 call cart39(flg1,flg2,gprimd,ipert2,natom,rprimd,vec1,vec2) 3188 do idir2=1,3 3189 d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2) 3190 ! And here carflg 3191 carflg(idir1,ipert1,idir2,ipert2)=flg2(idir2) 3192 end do 3193 end do 3194 end do 3195 end do 3196 end do 3197 3198 ! Second step 3199 do ipert1=1,mpert 3200 do ipert2=1,mpert 3201 do ii=1,2 3202 do idir2=1,3 3203 do idir1=1,3 3204 vec1(idir1)=d2cart(ii,idir1,ipert1,idir2,ipert2) 3205 ! Note here carflg 3206 flg1(idir1)=carflg(idir1,ipert1,idir2,ipert2) 3207 end do 3208 call cart39(flg1,flg2,gprimd,ipert1,natom,rprimd,vec1,vec2) 3209 do idir1=1,3 3210 d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1) 3211 ! And here carflg again 3212 carflg(idir1,ipert1,idir2,ipert2)=flg2(idir1) 3213 end do 3214 end do 3215 end do 3216 end do 3217 end do 3218 3219 end subroutine carteig2d
m_ddb/carttransf [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
carttransf
FUNCTION
Transform a second-derivative matrix (EIG2D) from reduced coordinates to cartesian coordinates.
INPUTS
blkflg(msize,nblok)=
( 1 if the element of the dynamical matrix has been calculated ;
0 otherwise )
gprimd(3,3)=basis vector in the reciprocal space
iqpt = number of the Q-point currently used
mband = maximal number of bands
mpert = maximum number of ipert
msize = size of the EIG2D arrays (3*mpert*3*mpert)
natom = number of atom
nblok = number of bloks in blkflg
nkpt = number of K-points
rprimd(3,3) = basis vector in the real space
OUTPUT
carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian EIG2D matrix has been calculated correctly ; 0 otherwise ) SIDE EFFECT blkval2(2,msize,mband,nkpt)=Second order eigenvalues (EIG2D) is transformed from reduced coordinates to cartesian coordinates
SOURCE
3053 subroutine carttransf(blkflg,blkval2,carflg,gprimd,iqpt,mband, mpert,msize,natom,nblok,nkpt,rprimd) 3054 3055 !Arguments ------------------------------- 3056 !scalars 3057 integer,intent(in) :: mband,msize 3058 integer,intent(in) :: iqpt 3059 integer,intent(in) :: mpert,nblok 3060 integer,intent(inout) :: natom,nkpt 3061 !arrays 3062 integer,intent(in) :: blkflg(msize,nblok) 3063 integer,intent(out) :: carflg(3,mpert,3,mpert) 3064 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 3065 real(dp),intent(inout) :: blkval2(2,msize,mband,nkpt) 3066 3067 !Local variables------------------------------- 3068 !scalars 3069 integer :: iatom1,iatom2,iband,idir1,idir2,ikpt, index 3070 !arrays 3071 real(dp),allocatable :: blkflgtmp(:,:,:,:,:), blkval2tmp(:,:,:,:,:,:), d2cart(:,:,:,:,:) 3072 3073 ! ********************************************************************* 3074 3075 ! Start by allocating local arrays 3076 ABI_MALLOC(blkflgtmp,(3,mpert,3,mpert,1)) 3077 ABI_MALLOC(blkval2tmp,(2,3,mpert,3,mpert,1)) 3078 ABI_MALLOC(d2cart,(2,3,mpert,3,mpert)) 3079 3080 ! Begin by formating the arrays to be compatible with cart29 3081 ! Then call cart29 to transform the arrays in cartesian coordinates 3082 ! Finally reformat the cartesian arrays in old format 3083 do ikpt=1,nkpt 3084 do iband=1,mband 3085 3086 do idir1=1,3 3087 do iatom1=1,mpert 3088 do idir2=1,3 3089 do iatom2=1,mpert 3090 index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1))) 3091 blkflgtmp(idir1,iatom1,idir2,iatom2,1) = blkflg(index,iqpt) 3092 blkval2tmp(:,idir1,iatom1,idir2,iatom2,1) = blkval2(:,index,iband,ikpt) 3093 end do 3094 end do 3095 end do 3096 end do 3097 3098 ! The 1sin the argument of cart29 are respectively iblok and nblok. We are doing only one blok. 3099 call carteig2d(blkflg(:,iqpt),blkval2tmp,carflg,d2cart,gprimd,1,mpert,natom,1,rprimd) 3100 3101 do idir1=1,3 3102 do iatom1=1,mpert 3103 do idir2=1,3 3104 do iatom2=1,mpert 3105 index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1))) 3106 blkval2(:,index,iband,ikpt) = d2cart(:,idir1,iatom1,idir2,iatom2) 3107 end do 3108 end do 3109 end do 3110 end do 3111 3112 end do 3113 end do 3114 3115 ABI_FREE(blkflgtmp) 3116 ABI_FREE(blkval2tmp) 3117 ABI_FREE(d2cart) 3118 3119 end subroutine carttransf
m_ddb/chkin9 [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
chkin9
FUNCTION
Check the value of some input parameters. Send error message and stop if needed. Also transform the meaning of atifc
INPUTS
atifc(natom)=list of the atom ifc to be analysed natifc= number of atom ifc to be analysed natom= number of atoms
OUTPUT
atifc(natom) = atifc(ia) equals 1 if the analysis of ifc has to be done for atom ia; otherwise 0.
NOTES
Only for one processor (no use of wrtout)
SOURCE
2185 subroutine chkin9(atifc,natifc,natom) 2186 2187 !Arguments ------------------------------- 2188 !scalars 2189 integer,intent(in) :: natifc,natom 2190 !arrays 2191 integer,intent(inout) :: atifc(natom) 2192 2193 !Local variables ------------------------- 2194 !scalars 2195 integer :: iatifc 2196 character(len=500) :: msg 2197 !arrays 2198 integer,allocatable :: work(:) 2199 2200 ! ********************************************************************* 2201 2202 if(natifc>natom)then 2203 write(msg, '(a,i0,3a,i0,3a)' )& 2204 'The number of atom ifc in the input files',natifc,',',ch10,& 2205 'is larger than the number of atoms',natom,'.',ch10,& 2206 'Action: change natifc in the input file.' 2207 ABI_ERROR(msg) 2208 end if 2209 2210 if(natifc>=1)then 2211 ABI_MALLOC(work,(natom)) 2212 work(:)=0 2213 2214 do iatifc=1,natifc 2215 if(atifc(iatifc)<=0.or.atifc(iatifc)>natom)then 2216 write(msg, '(a,i0,5a,i0,3a)' )& 2217 'For iatifc=',iatifc,', the number of the atom ifc to be ',ch10,& 2218 'analysed is not valid : either negative, ',ch10,& 2219 'zero, or larger than natom =',natom,'.',ch10,& 2220 'Action: change atifc in your input file.' 2221 ABI_ERROR(msg) 2222 end if 2223 work(atifc(iatifc))=1 2224 end do 2225 2226 atifc(1:natom)=work(:) 2227 ABI_FREE(work) 2228 end if 2229 2230 end subroutine chkin9
m_ddb/ddb_bcast [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_bcast
FUNCTION
MPI broadcast all types for the ddb_type structure
INPUTS
comm=MPI communicator
SIDE EFFECTS
Ddb<type(ddb_type)>= Input if node is master, other nodes returns with a completely initialized instance.
SOURCE
1249 subroutine ddb_bcast(ddb, comm) 1250 1251 !Arguments ------------------------------- 1252 !array 1253 class(ddb_type),intent(inout) :: ddb 1254 integer, intent(in) :: comm 1255 1256 !Local variables------------------------------- 1257 !scalars 1258 integer, parameter :: master=0 1259 integer :: ierr 1260 ! ************************************************************************* 1261 1262 if (xmpi_comm_size(comm) == 1) return 1263 1264 DBG_ENTER("COLL") 1265 1266 ! Transmit dimensions and static variables. 1267 call xmpi_bcast(ddb%nblok, master, comm, ierr) 1268 call xmpi_bcast(ddb%natom, master, comm, ierr) 1269 call xmpi_bcast(ddb%ntypat, master, comm, ierr) 1270 call xmpi_bcast(ddb%nsppol, master, comm, ierr) 1271 call xmpi_bcast(ddb%mpert, master, comm, ierr) 1272 call xmpi_bcast(ddb%msize, master, comm, ierr) 1273 1274 call xmpi_bcast(ddb%occopt, master, comm, ierr) 1275 call xmpi_bcast(ddb%prtvol, master, comm, ierr) 1276 1277 !real 1278 call xmpi_bcast(ddb%rprim, master, comm, ierr) 1279 call xmpi_bcast(ddb%gprim, master, comm, ierr) 1280 call xmpi_bcast(ddb%acell, master, comm, ierr) 1281 1282 ! Allocate arrays on the other nodes. 1283 if (xmpi_comm_rank(comm) /= master) then 1284 call ddb%malloc(ddb%msize, ddb%nblok, ddb%natom, ddb%ntypat, ddb%mpert) 1285 end if 1286 1287 call xmpi_bcast(ddb%flg, master, comm, ierr) 1288 call xmpi_bcast(ddb%typ, master, comm, ierr) 1289 call xmpi_bcast(ddb%amu, master, comm, ierr) 1290 call xmpi_bcast(ddb%nrm, master, comm, ierr) 1291 call xmpi_bcast(ddb%qpt, master, comm, ierr) 1292 call xmpi_bcast(ddb%val, master, comm, ierr) 1293 1294 DBG_EXIT("COLL") 1295 1296 end subroutine ddb_bcast
m_ddb/ddb_can_merge_blocks [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_can_merge_blocks
FUNCTION
Return true if iblok1 of ddb1 can be merged to iblok2 of ddb2
INPUTS
ddb1=ddb object 1 ddb2=ddb object 2 iblok1=block index from ddb1 iblok2=block index from ddb2
OUTPUT
can_merge=.true. if the blocks are compatible for merging.
SOURCE
2801 logical function ddb_can_merge_blocks(ddb1, ddb2, iblok1, iblok2) result(can_merge) 2802 2803 !Arguments ------------------------------- 2804 !array 2805 class(ddb_type),intent(inout) :: ddb1 2806 type(ddb_type),intent(inout) :: ddb2 2807 integer,intent(in) :: iblok1 2808 integer,intent(in) :: iblok2 2809 2810 !local variables 2811 !scalars 2812 integer :: nq, ii, blktyp 2813 real(dp),parameter :: qtol=2.0d-8 2814 real(dp) :: diff 2815 2816 ! ************************************************************************ 2817 2818 can_merge = .false. 2819 if(ddb1%typ(iblok1)/=ddb2%typ(iblok2)) return 2820 2821 blktyp = ddb1%typ(iblok1) 2822 2823 can_merge = .true. 2824 2825 if (is_type_d0E(blktyp) .or. is_type_d1E(blktyp)) return 2826 2827 ! Compare wavevectors 2828 if (is_type_d2E(blktyp).or.is_type_d2eig(blktyp))then 2829 nq=1 2830 else if (is_type_d3E(blktyp))then 2831 nq=3 2832 end if 2833 2834 do ii=1,nq 2835 diff = (ddb1%qpt(1+3*(ii-1),iblok1)/ddb1%nrm(ii,iblok1) & 2836 - ddb2%qpt(1+3*(ii-1),iblok2)/ddb2%nrm(ii,iblok2)) 2837 if (abs(diff) > qtol) can_merge = .false. 2838 diff = (ddb1%qpt(2+3*(ii-1),iblok1)/ddb1%nrm(ii,iblok1) & 2839 - ddb2%qpt(2+3*(ii-1),iblok2)/ddb2%nrm(ii,iblok2)) 2840 if (abs(diff) > qtol) can_merge = .false. 2841 diff = (ddb1%qpt(3+3*(ii-1),iblok1)/ddb1%nrm(ii,iblok1) & 2842 - ddb2%qpt(3+3*(ii-1),iblok2)/ddb2%nrm(ii,iblok2)) 2843 if (abs(diff) > qtol) can_merge = .false. 2844 end do 2845 2846 end function ddb_can_merge_blocks
m_ddb/ddb_copy [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_copy
FUNCTION
Create object and copy all types for the ddb_type structure
SOURCE
578 subroutine ddb_copy(iddb, oddb) 579 580 !Arguments ------------------------------- 581 !array 582 class(ddb_type),intent(in) :: iddb 583 type(ddb_type),intent(out) :: oddb 584 585 ! ************************************************************************ 586 587 ! Copy dimensions and static variables. 588 oddb%msize = iddb%msize 589 oddb%mpert = iddb%mpert 590 oddb%nblok = iddb%nblok 591 oddb%natom = iddb%natom 592 oddb%ntypat = iddb%ntypat 593 oddb%occopt = iddb%occopt 594 oddb%prtvol = iddb%prtvol 595 596 oddb%rprim = iddb%rprim 597 oddb%gprim = iddb%gprim 598 oddb%acell = iddb%acell 599 600 ! Allocate and copy the allocatable arrays. 601 call alloc_copy(iddb%flg, oddb%flg) 602 call alloc_copy(iddb%typ, oddb%typ) 603 call alloc_copy(iddb%amu, oddb%amu) 604 call alloc_copy(iddb%nrm, oddb%nrm) 605 call alloc_copy(iddb%qpt, oddb%qpt) 606 call alloc_copy(iddb%val, oddb%val) 607 608 end subroutine ddb_copy
m_ddb/ddb_diagoq [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_diagoq
FUNCTION
Compute the phonon frequencies at the specified q-point by performing a direct diagonalization of the dynamical matrix. The q-point **MUST** be one the points stored in the DDB file.
INPUTS
ddb<type(ddb_type)>=Object storing the DDB results. crystal<type(crystal_t)> = Information on the crystalline structure. asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file. symdynmat=If equal to 1, the dynamical matrix is symmetrized in dfpt_phfrq before the diagonalization. rftyp = 1 if non-stationary block 2 if stationary block 3 if third order derivatives qpt(3)=q-point in reduced coordinates.
OUTPUT
phfrq(3*crystal%natom)=Phonon frequencies in Hartree displ_cart(2,3*%natom,3*%natom)=Phonon displacement in Cartesian coordinates [out_eigvec(2*3*natom*3*natom) = The igenvectors of the dynamical matrix. [out_displ_red(2*3*natom*3*natom) = The displacement in reduced coordinates.
SOURCE
4399 subroutine ddb_diagoq(ddb, crystal, qpt, asrq0, symdynmat, rftyp, phfrq, displ_cart, & 4400 out_eigvec,out_displ_red) ! Optional [out] 4401 4402 !Arguments ------------------------------- 4403 !scalars 4404 integer,intent(in) :: symdynmat 4405 integer,intent(in) :: rftyp 4406 class(ddb_type),intent(in) :: ddb 4407 type(asrq0_t),intent(inout) :: asrq0 4408 type(crystal_t),intent(in) :: crystal 4409 !arrays 4410 real(dp),intent(in) :: qpt(3) 4411 real(dp),intent(out) :: displ_cart(2,3,crystal%natom,3,crystal%natom) 4412 real(dp),intent(out) :: phfrq(3*crystal%natom) 4413 !real(dp),optional,intent(out) :: out_d2cart(2,3*crystal%natom,3*crystal%natom) 4414 real(dp),optional,intent(out) :: out_eigvec(2,3,crystal%natom,3*crystal%natom) 4415 real(dp),optional,intent(out) :: out_displ_red(2,3,crystal%natom,3*crystal%natom) 4416 4417 !Local variables------------------------------- 4418 integer :: iblok,natom 4419 !arrays 4420 integer :: rfphon(4),rfelfd(4),rfstrs(4) 4421 real(dp) :: qphnrm(3), qphon_padded(3,3),d2cart(2,ddb%msize),my_qpt(3) 4422 real(dp) :: eigvec(2,3,crystal%natom,3*crystal%natom),eigval(3*crystal%natom) 4423 4424 ! ************************************************************************ 4425 4426 ! Use my_qpt because dfpt_phfrq can change the q-point (very bad design) 4427 qphnrm = one; my_qpt = qpt 4428 4429 ! Look for the information in the DDB (no interpolation here!) 4430 rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0 4431 qphon_padded = zero; qphon_padded(:,1) = qpt 4432 natom = crystal%natom 4433 4434 call ddb%get_block(iblok, qphon_padded, qphnrm, rfphon, rfelfd, rfstrs, rftyp) 4435 ABI_CHECK(iblok /= 0, sjoin("Cannot find q-point ", ktoa(qpt)," in DDB file")) 4436 4437 ! Copy the dynamical matrix in d2cart 4438 d2cart(:,1:ddb%msize) = ddb%val(:,:,iblok) 4439 4440 ! Eventually impose the acoustic sum rule based on previously calculated d2asr 4441 call asrq0%apply(natom, ddb%mpert, ddb%msize, crystal%xcart, d2cart) 4442 4443 ! Calculation of the eigenvectors and eigenvalues of the dynamical matrix 4444 call dfpt_phfrq(ddb%amu,displ_cart,d2cart,eigval,eigvec,crystal%indsym,& 4445 ddb%mpert,crystal%nsym,natom,crystal%nsym,crystal%ntypat,phfrq,qphnrm(1),my_qpt,& 4446 crystal%rprimd,symdynmat,crystal%symrel,crystal%symafm,crystal%typat,crystal%ucvol) 4447 4448 ! Return the dynamical matrix and the eigenvector for this q-point 4449 !if (present(out_d2cart)) out_d2cart = d2cart(:,:3*natom,:3*natom) 4450 if (present(out_eigvec)) out_eigvec = eigvec 4451 4452 ! Return phonon displacement in reduced coordinates. 4453 if (present(out_displ_red)) call phdispl_cart2red(natom, crystal%gprimd, displ_cart, out_displ_red) 4454 4455 end subroutine ddb_diagoq
m_ddb/ddb_free [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_free
FUNCTION
Clean and deallocate types for the ddb_type structure
SOURCE
545 subroutine ddb_free(ddb) 546 547 !Arguments ------------------------------- 548 class(ddb_type),intent(inout) :: ddb 549 550 ! ************************************************************************ 551 552 !integer 553 ABI_SFREE(ddb%flg) 554 ABI_SFREE(ddb%typ) 555 556 ! real 557 ABI_SFREE(ddb%amu) 558 ABI_SFREE(ddb%qpt) 559 ABI_SFREE(ddb%nrm) 560 ABI_SFREE(ddb%kpt) 561 ABI_SFREE(ddb%val) 562 ABI_SFREE(ddb%eig2dval) 563 564 end subroutine ddb_free
m_ddb/ddb_from_file [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_from_file
FUNCTION
This subroutine reads data from the DDB file and constructs an instance of ddb_type It also returns an instance of crystal_t with the crystalline structure reported in the DDB file and the DDB header.
INPUTS
filename=DDB filename. comm=MPI communicator. [prtvol] = Verbosity level raw = 1 -> do not perform any symetrization or transformation to cartesian coordinates. 0 (default) -> do perform these transformations.
OUTPUT
ddb<type(ddb_type)>=Object storing the DDB results. crystal<type(crystal_t)>=Crystal structure parameters ddb_hdr<type(ddb_hdr_type)>= Header of the DDB file.
SOURCE
2399 subroutine ddb_from_file(ddb, filename, ddb_hdr, crystal, comm, prtvol, raw) 2400 2401 !Arguments ------------------------------- 2402 !scalars 2403 class(ddb_type),intent(inout) :: ddb 2404 integer,intent(in) :: comm 2405 integer,optional,intent(in) :: prtvol, raw 2406 character(len=*),intent(in) :: filename 2407 type(crystal_t),intent(out) :: Crystal 2408 type(ddb_hdr_type),intent(out) :: ddb_hdr 2409 !array 2410 2411 !Local variables------------------------------- 2412 integer :: iomode 2413 character(len=fnlen) :: filename_ 2414 integer :: prtvol_ 2415 character(len=500) :: msg 2416 2417 ! ************************************************************************ 2418 2419 DBG_ENTER("COLL") 2420 2421 prtvol_ = 0; if (present(prtvol)) prtvol_ = prtvol 2422 2423 call ddb_hdr%get_iomode(filename, 1, iomode, filename_) 2424 2425 if (iomode==IO_MODE_ETSF) then 2426 call ddb%read_nc(filename_, ddb_hdr, crystal, comm, prtvol, raw) 2427 else if (iomode==IO_MODE_FORTRAN) then 2428 call ddb%read_txt(filename_, ddb_hdr, crystal, comm, prtvol, raw) 2429 end if 2430 2431 ! Print out info on the crystal 2432 if (prtvol_ >= 0) then 2433 2434 call ddb_hdr%crystal%print(unit=ab_out) 2435 call ddb_hdr%crystal%print(unit=std_out) 2436 2437 write(msg, '(2a,i0,a)' )ch10,' DDB file with ',ddb%nblok,' blocks has been read.' 2438 call wrtout(std_out,msg) 2439 call wrtout(ab_out,msg) 2440 2441 end if 2442 2443 DBG_EXIT("COLL") 2444 2445 end subroutine ddb_from_file
m_ddb/ddb_get_asrq0 [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_asrq0
FUNCTION
In case the interatomic forces are not calculated, the ASR-correction has to be determined here from the Dynamical matrix at Gamma. In case the DDB does not contain this information, the subroutine returns iblok=0 %d2asr is initialized and set to zero to preserve the old behaviour.
INPUTS
asr=Input variable selecting the method for the ASR
rftyp = 1 if non-stationary block
2 if stationary block
3 if third order derivatives
xcart(3,ddb%atom)=Cartesian coordinates of the atoms.
SIDE EFFECTS
ddb<type(ddb_type)>= Database with the derivates. The routine does not change it except when asr is in [3,4]. TODO This should not happen.
OUTPUT
asrq0<asrq0_t> iblok= is set to 0 if the Gamma block is not found
SOURCE
4106 type(asrq0_t) function ddb_get_asrq0(ddb, asr, rftyp, xcart) result(asrq0) 4107 4108 !Arguments ------------------------------- 4109 !scalars 4110 integer,intent(in) :: asr,rftyp 4111 class(ddb_type),intent(inout) :: ddb 4112 !arrays 4113 real(dp),intent(in) :: xcart(3,ddb%natom) 4114 4115 !Local variables------------------------------- 4116 !scalars 4117 integer :: dims,iblok 4118 !character(len=500) :: msg 4119 !arrays 4120 integer :: rfelfd(4),rfphon(4),rfstrs(4) 4121 real(dp) :: qphnrm(3),qphon(3,3) 4122 real(dp),allocatable :: d2asr_res(:,:,:,:,:),d2cart(:,:) 4123 4124 ! ************************************************************************ 4125 4126 asrq0%asr = asr; asrq0%natom = ddb%natom 4127 4128 ! Find the Gamma block in the DDB (no need for E-field entries) 4129 qphon(:,1)=zero 4130 qphnrm(1)=zero 4131 rfphon(1:2)=1 4132 rfelfd(:)=0 4133 rfstrs(:)=0 4134 4135 call ddb%get_block(asrq0%iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 4136 ! this is to maintain the old behaviour in which the arrays where allocated and set to zero in anaddb. 4137 ABI_MALLOC(asrq0%d2asr, (2,3,ddb%natom,3,ddb%natom)) 4138 asrq0%d2asr = zero 4139 4140 if (asrq0%iblok == 0) return 4141 iblok = asrq0%iblok 4142 4143 select case (asrq0%asr) 4144 case (0) 4145 continue 4146 4147 case (1,2) 4148 call asria_calc(asr,asrq0%d2asr,ddb%val(:,:,iblok),ddb%mpert,ddb%natom) 4149 4150 case (3,4) 4151 ! Rotational invariance for 1D and 0D systems 4152 ! Compute uinvers, vtinvers and singular matrices. 4153 dims = 3*ddb%natom*(3*ddb%natom-1) / 2 4154 ABI_CALLOC(asrq0%uinvers, (dims, dims)) 4155 ABI_CALLOC(asrq0%vtinvers,(dims, dims)) 4156 ABI_CALLOC(asrq0%singular, (dims)) 4157 4158 call asrprs(asr,1,3,asrq0%uinvers,asrq0%vtinvers,asrq0%singular,& 4159 ddb%val(:,:,iblok),ddb%mpert,ddb%natom,xcart) 4160 4161 case (5) 4162 ! d2cart is a temp variable here 4163 ABI_MALLOC(d2cart,(2,ddb%msize)) 4164 d2cart = ddb%val(:,:,iblok) 4165 ! calculate diagonal correction 4166 call asria_calc(2,asrq0%d2asr,d2cart,ddb%mpert,ddb%natom) 4167 ! apply diagonal correction 4168 call asria_corr(2,asrq0%d2asr,d2cart,ddb%mpert,ddb%natom) 4169 ! hermitianize 4170 call mkherm(d2cart,3*ddb%mpert) 4171 ! remove remaining ASR rupture due to Hermitianization 4172 ABI_MALLOC(d2asr_res,(2,3,ddb%natom,3,ddb%natom)) 4173 call asria_calc(asr,d2asr_res,d2cart,ddb%mpert,ddb%natom) 4174 ! full correction is sum of both 4175 asrq0%d2asr = asrq0%d2asr + d2asr_res 4176 4177 ABI_FREE(d2cart) 4178 ABI_FREE(d2asr_res) 4179 4180 case default 4181 ABI_ERROR(sjoin("Wrong value for asr:", itoa(asr))) 4182 end select 4183 4184 end function ddb_get_asrq0
m_ddb/ddb_get_block [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_block
FUNCTION
This routine finds the block that contains the information on the derivatives of the total energy specified by the parameters rfphon,rfelfd,rfstrs,rftyp and the phonon wavevectors qphon (and their normalisation). In case the DDB does not contain this information, the subroutine returns iblok=0
INPUTS
ddb = ddb blok datastructure
flg(msize,nblok)=flag for every matrix element:
0 => the element is not in the data block.
1 => the element is in the data blok.
nrm(3,nblok)=normalization factors for the three allowed wavevectors
qpt(3,nblok)=wavevector of the perturbation(s). The elements
typ(nblok)=type of the block.
(1=> non-stationary block),
(2=> stationary block),
(3=> third order derivative).
qphon(3,3)=wavevectors for the three possible phonons
(note : only one should be used in case of second derivative of total energy,
because we know that the second is the opposite of this value)
qphnrm(3) =normalisation factors for the three possible phonons
rfphon(4) = 1=> response to phonons
2=> second derivative of total energy
rfelfd(4) = 1=> d/dk, 2=> electric field only, 3=> both (see comment on rfphon)
rfstrs(4) = 1=> uniaxial stresses, 2=> shear stresses, 3=> both (see comment on rfphon)
rftyp =
0 => total energy
1 => non-stationary formulation of the 2nd derivative
2 => stationary formulation of the 2nd derivative
3 => third derivative of total energy
4 => first-order derivatives of total energy
33 => long wave third order derivatives of total energy
85 => molecular Berry curvature
[rfqvec(4)] = 1=> d/dq (optional)
OUTPUT
iblok= number of the block that corresponds to the specifications. 0 if not found.
SOURCE
1347 subroutine ddb_get_block(ddb, iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp, rfqvec) 1348 1349 !Arguments ------------------------------- 1350 !scalars 1351 integer,intent(in) :: rftyp 1352 integer,intent(out) :: iblok 1353 class(ddb_type),intent(in) :: ddb 1354 !arrays 1355 integer,intent(in) :: rfelfd(4),rfphon(4),rfstrs(4) 1356 real(dp),intent(inout) :: qphnrm(3),qphon(3,3) 1357 integer,optional,intent(in) :: rfqvec(4) 1358 1359 !Local variables ------------------------- 1360 !scalars 1361 integer :: blkgam,ider,idir,idir1,idir2,idir3,ii,index,ipert,ipert1,ipert2 1362 integer :: ipert3,nder,ok,mpert,natom 1363 character(len=500) :: msg 1364 !arrays 1365 integer :: gamma(3) 1366 integer,allocatable :: worki(:,:) 1367 real(dp) :: qpt(3) 1368 integer :: rfqvec_(4) 1369 1370 ! ********************************************************************* 1371 1372 mpert = ddb%mpert 1373 natom = ddb%natom 1374 1375 ! Get the number of derivative 1376 if (is_type_d2E(rftyp)) then 1377 nder=2 1378 else if (is_type_d3E(rftyp)) then 1379 nder=3 1380 else if (is_type_d0E(rftyp)) then 1381 nder=0 1382 else if (is_type_d1E(rftyp)) then 1383 nder=1 1384 else 1385 write(msg, '(a,i0,a)')' rftyp is equal to ',rftyp,'. The only allowed values are 0, 1, 2, 3, 5, 6, 33 or 85.' 1386 ABI_BUG(msg) 1387 end if 1388 1389 rfqvec_(:)=0; if(present(rfqvec))rfqvec_(:)=rfqvec(:) 1390 1391 ! In case of a second-derivative, a second phonon wavevector is provided. 1392 if(nder==2)then 1393 do ii=1,3 1394 qphon(ii,2)=-qphon(ii,1) 1395 end do 1396 qphnrm(2)=qphnrm(1) 1397 end if 1398 1399 ! In case of a third derivative, the sum of wavevectors to gamma is checked 1400 if (nder == 3) then 1401 qpt(:) = qphon(:,1)/qphnrm(1) + qphon(:,2)/qphnrm(2) + qphon(:,3)/qphnrm(3) 1402 call gamma9(gamma(nder),qpt,qphnrm(1),DDB_QTOL) 1403 if (gamma(nder) == 0) then 1404 write(msg,'(a,a,a)')& 1405 'the sum of the wavevectors of the third-order energy is ',ch10,& 1406 'not equal to zero' 1407 ABI_ERROR(msg) 1408 end if 1409 end if 1410 1411 ! Check the validity of the requirement 1412 do ider=1,nder 1413 ! Identifies if qphon is at gamma 1414 call gamma9(gamma(ider),qphon(1:3,ider),qphnrm(ider),DDB_QTOL) 1415 1416 if(gamma(ider)==0)then 1417 if(rfstrs(ider)/=0.or.rfelfd(ider)/=0.or.rfqvec_(ider)/=0)then 1418 write(msg, '(a,a)' )& 1419 'Not yet able to handle stresses or electric fields',ch10,& 1420 'with non-zero wavevector.' 1421 ABI_BUG(msg) 1422 end if 1423 end if 1424 end do 1425 1426 ! Initialise the perturbation table 1427 ABI_MALLOC(worki,(mpert,4)) 1428 worki(:,1:nder)=0 1429 1430 ! Build the perturbation table 1431 do ider=1,nder 1432 ! First the phonons 1433 if(rfphon(ider)==1)then ! what about the rfphon = 2 case ? 1434 do ipert=1,natom 1435 worki(ipert,ider)=1 1436 end do 1437 end if 1438 ! Then the d/dk 1439 if (rfelfd(ider)==1.or.rfelfd(ider)==3) worki(natom+1,ider)=1 1440 ! Then the electric field 1441 if (rfelfd(ider)==2.or.rfelfd(ider)==3) worki(natom+2,ider)=1 1442 ! Then the ddq 1443 if (rfqvec_(ider)==1) worki(natom+8,ider)=1 1444 ! Then the uniaxial stress 1445 if (rfstrs(ider)==1.or.rfstrs(ider)==3) worki(natom+3,ider)=1 1446 ! At last, the shear stress 1447 if(rfstrs(ider)==2.or.rfstrs(ider)==3) worki(natom+4,ider)=1 1448 end do 1449 1450 ! Examine every blok: 1451 do iblok=1,ddb%nblok 1452 1453 ! If this variable is still 1 at the end of the examination, the blok is the good one... 1454 ok=1 1455 1456 ! Check the type 1457 if(rftyp/=ddb%typ(iblok)) ok=0 1458 1459 ! Check the wavevector 1460 if( ok==1 )then 1461 1462 if (nder == 2) then 1463 call gamma9(blkgam,ddb%qpt(1:3,iblok),ddb%nrm(1,iblok),DDB_QTOL) 1464 if(blkgam/=gamma(1))then 1465 ok=0 1466 else if(blkgam==0)then 1467 do idir=1,3 1468 if( abs( ddb%qpt(idir,iblok)/ddb%nrm(1,iblok) - qphon(idir,1)/qphnrm(1) )>DDB_QTOL ) ok=0 1469 end do 1470 end if 1471 1472 else if (nder == 3) then 1473 do ider = 1, nder 1474 do idir=1,3 1475 if( abs( ddb%qpt(idir+3*(ider-1),iblok)/ddb%nrm(ider,iblok) - qphon(idir,ider)/qphnrm(ider) )>DDB_QTOL )then 1476 ok=0 1477 end if ! qphon 1478 end do ! idir 1479 end do ! nder 1480 end if ! nder 1481 1482 end if ! ok 1483 1484 ! Check if there is enough information in this blok 1485 if( ok==1 )then 1486 1487 if (nder == 0) then 1488 if (ddb%flg(1,iblok) /= 1) then 1489 ok = 0 1490 if (ddb%prtvol > 1) then 1491 write(msg,'(a,i0,3a)' )& 1492 'The block ',iblok,' does not match the requirement',ch10,& 1493 'because it lacks the total energy' 1494 ABI_COMMENT(msg) 1495 end if 1496 end if 1497 end if 1498 1499 do ipert1=1,mpert 1500 1501 if ((nder == 4).and.(worki(ipert1,4) == 1).and.(ok == 1)) then 1502 do idir1 = 1, 3 1503 index = 3*(ipert1 - 1) + idir1 1504 if (ddb%flg(index,iblok) /= 1) ok = 0 1505 end do 1506 end if 1507 1508 if (worki(ipert1,1)==1 .and. ok==1 )then 1509 do ipert2=1,mpert 1510 if (worki(ipert2,2)==1 .and. ok==1 )then 1511 do idir1=1,3 1512 do idir2=1,3 1513 1514 if (nder == 2) then 1515 index=idir1+ 3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1))) 1516 if (ddb%flg(index,iblok)/=1) ok=0 1517 1518 else if (nder == 3) then 1519 do ipert3 = 1, mpert 1520 if (worki(ipert3,3) == 1 .and. ok == 1) then 1521 do idir3 = 1, 3 1522 index = idir1 + & 1523 3*((ipert1 - 1) + mpert*((idir2 - 1) + & 1524 3*((ipert2 -1 ) + mpert*((idir3 - 1) + 3*(ipert3 - 1))))) 1525 if (ddb%flg(index,iblok) /= 1) ok = 0 1526 end do ! idir3 1527 end if ! worki(ipert3,3) 1528 end do ! i3pert 1529 end if 1530 1531 end do 1532 end do 1533 end if 1534 end do 1535 end if 1536 end do 1537 end if 1538 1539 ! Now that everything has been checked, eventually end the search 1540 if(ok==1)exit 1541 end do 1542 1543 if(ok==0)then 1544 iblok=0 1545 1546 if (ddb%prtvol > 1) then 1547 write(msg, '(3a)' )& 1548 ' gtblk9 : ',ch10,& 1549 ' Unable to find block corresponding to the following specifications :' 1550 call wrtout(std_out,msg) 1551 write(msg, '(a,i3)' )' Type (rfmeth) =',rftyp 1552 call wrtout(std_out,msg) 1553 write(msg, '(a)' ) ' ider qphon(3) qphnrm rfphon rfelfd rfstrs rfqvec' 1554 call wrtout(std_out,msg) 1555 do ider=1,nder 1556 write(msg, '(i4,4f6.2,4i7)' )& 1557 ider,(qphon(ii,ider),ii=1,3),qphnrm(ider),rfphon(ider),rfelfd(ider),rfstrs(ider),rfqvec_(ider) 1558 call wrtout(std_out,msg) 1559 end do 1560 end if 1561 end if 1562 1563 if (ok==1 .and. ddb%prtvol > 1) then 1564 write(msg,'(a,i0,2a)')' gtblk9: found block number ',iblok,' agree with',' specifications ' 1565 call wrtout(std_out,msg) 1566 end if 1567 1568 ABI_FREE(worki) 1569 1570 end subroutine ddb_get_block
m_ddb/ddb_get_d1matr [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_d1matr
FUNCTION
Transform the first-order derivative matrix from flat indices to real tensor d1matr(cplex,ncart,natom)
INPUTS
iblok=index of the block to get.
OUTPUT
d1matr=the first-order derivative matrix. flg=flag to indicate presence of a given element.
SOURCE
1023 subroutine ddb_get_d1matr(ddb, iblok, d1matr, flg) 1024 1025 !Arguments ------------------------------- 1026 !array 1027 class(ddb_type),intent(inout) :: ddb 1028 integer,intent(in) :: iblok 1029 real(dp), allocatable, intent(out) :: d1matr(:,:,:) 1030 integer, allocatable, intent(out) :: flg(:,:) 1031 !scalars 1032 1033 !Local variables ------------------------- 1034 !scalars 1035 integer :: ii,idir1,ipert1 1036 1037 ! ************************************************************************ 1038 1039 ABI_MALLOC(d1matr, (2,3,ddb%mpert)) 1040 ABI_MALLOC(flg, (3,ddb%mpert)) 1041 1042 d1matr = zero 1043 1044 ii=0 1045 do ipert1=1,ddb%mpert 1046 do idir1=1,3 1047 ii=ii+1 1048 flg(idir1,ipert1) = ddb%flg(ii,iblok) 1049 if (ddb%flg(ii,iblok) > 0) then 1050 d1matr(1,idir1,ipert1) = ddb%val(1,ii,iblok) 1051 d1matr(2,idir1,ipert1) = ddb%val(2,ii,iblok) 1052 end if 1053 end do 1054 end do 1055 1056 end subroutine ddb_get_d1matr
m_ddb/ddb_get_d2eig [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_d2eig
FUNCTION
Transform the second-order derivative matrix of eigenvalues from flat indices to real tensor d2eig(cplex,ncart,natom,ncart,natom,nband,nkpt)
INPUTS
iblok=index of the block to get.
OUTPUTS
d2eig(cplex,ncart,natom,ncart,natom,nband,nkpt) =
the second-order derivative matrix of eigenvalues
with the isppol index wrapped in nband.
flg(ncart,natom,ncart,natom)=flag to indicate presence of a given element.
SOURCE
5820 subroutine ddb_get_d2eig(ddb, d2eig, flg, iblok) 5821 5822 !Arguments ------------------------------- 5823 !array 5824 class(ddb_type),intent(inout) :: ddb 5825 real(dp), allocatable, intent(out) :: d2eig(:,:,:,:,:,:,:) 5826 integer,intent(in) :: iblok 5827 integer, allocatable, intent(out) :: flg(:,:,:,:) 5828 !scalars 5829 5830 !Local variables ------------------------- 5831 !scalars 5832 integer :: ii,idir1,idir2,ipert1,ipert2,iband,ikpt 5833 5834 ! ************************************************************************ 5835 5836 ABI_MALLOC(d2eig, (2,3,ddb%mpert,3,ddb%mpert,ddb%nband,ddb%nkpt)) 5837 ABI_MALLOC(flg, (3,ddb%mpert,3,ddb%mpert)) 5838 5839 d2eig = zero 5840 5841 do ikpt=1,ddb%nkpt 5842 do iband=1,ddb%nband 5843 ii=0 5844 do ipert2=1,ddb%mpert 5845 do idir2=1,3 5846 do ipert1=1,ddb%mpert 5847 do idir1=1,3 5848 ii=ii+1 5849 flg(idir1,ipert1,idir2,ipert2) = ddb%flg(ii,iblok) 5850 if (ddb%flg(ii,iblok) > 0) then 5851 d2eig(1,idir1,ipert1,idir2,ipert2,iband,ikpt) = ddb%eig2dval(1,ii,iband,ikpt) 5852 d2eig(2,idir1,ipert1,idir2,ipert2,iband,ikpt) = ddb%eig2dval(2,ii,iband,ikpt) 5853 end if 5854 end do 5855 end do 5856 end do 5857 end do 5858 end do 5859 end do 5860 5861 end subroutine ddb_get_d2eig
m_ddb/ddb_get_d2matr [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_d2matr
FUNCTION
Transform the second-order derivative matrix from flat indices to real tensor d2matr(cplex,ncart,natom,ncart,natom)
INPUTS
iblok=index of the block to get.
OUTPUT
d2matr=the second-order derivative matrix. flg=flag to indicate presence of a given element.
SOURCE
829 subroutine ddb_get_d2matr(ddb, iblok, d2matr, flg) 830 831 !Arguments ------------------------------- 832 !array 833 class(ddb_type),intent(inout) :: ddb 834 integer,intent(in) :: iblok 835 real(dp), allocatable, intent(out) :: d2matr(:,:,:,:,:) 836 integer, allocatable, intent(out) :: flg(:,:,:,:) 837 !scalars 838 839 !Local variables ------------------------- 840 !scalars 841 integer :: ii,idir1,idir2,ipert1,ipert2 842 843 ! ************************************************************************ 844 845 ABI_MALLOC(d2matr, (2,3,ddb%mpert,3,ddb%mpert)) 846 ABI_MALLOC(flg, (3,ddb%mpert,3,ddb%mpert)) 847 848 d2matr = zero 849 850 ii=0 851 do ipert2=1,ddb%mpert 852 do idir2=1,3 853 do ipert1=1,ddb%mpert 854 do idir1=1,3 855 ii=ii+1 856 flg(idir1,ipert1,idir2,ipert2) = ddb%flg(ii,iblok) 857 if (ddb%flg(ii,iblok) > 0) then 858 d2matr(1,idir1,ipert1,idir2,ipert2) = ddb%val(1,ii,iblok) 859 d2matr(2,idir1,ipert1,idir2,ipert2) = ddb%val(2,ii,iblok) 860 end if 861 end do 862 end do 863 end do 864 end do 865 866 end subroutine ddb_get_d2matr
m_ddb/ddb_get_d3matr [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_d3matr
FUNCTION
Transform the third-order derivative matrix from flat indices to real tensor d3matr(cplex,ncart,natom,ncart,natom,ncart,natom)
INPUTS
iblok=index of the block to get.
OUTPUT
d3matr=the third-order derivative matrix. flg=flag to indicate presence of a given element.
SOURCE
5755 subroutine ddb_get_d3matr(ddb, iblok, d3matr, flg) 5756 5757 !Arguments ------------------------------- 5758 !array 5759 class(ddb_type),intent(inout) :: ddb 5760 integer,intent(in) :: iblok 5761 real(dp), allocatable, intent(out) :: d3matr(:,:,:,:,:,:,:) 5762 integer, allocatable, intent(out) :: flg(:,:,:,:,:,:) 5763 !scalars 5764 5765 !Local variables ------------------------- 5766 !scalars 5767 integer :: ii,idir1,idir2,idir3,ipert1,ipert2,ipert3 5768 5769 ! ************************************************************************ 5770 5771 ABI_MALLOC(d3matr, (2,3,ddb%mpert,3,ddb%mpert,3,ddb%mpert)) 5772 ABI_MALLOC(flg, (3,ddb%mpert,3,ddb%mpert,3,ddb%mpert)) 5773 5774 d3matr = zero 5775 5776 ii=0 5777 do ipert3=1,ddb%mpert 5778 do idir3=1,3 5779 do ipert2=1,ddb%mpert 5780 do idir2=1,3 5781 do ipert1=1,ddb%mpert 5782 do idir1=1,3 5783 ii=ii+1 5784 flg(idir1,ipert1,idir2,ipert2,idir3,ipert3) = ddb%flg(ii,iblok) 5785 if (ddb%flg(ii,iblok) > 0) then 5786 d3matr(1,idir1,ipert1,idir2,ipert2,idir3,ipert3) = ddb%val(1,ii,iblok) 5787 d3matr(2,idir1,ipert1,idir2,ipert2,idir3,ipert3) = ddb%val(2,ii,iblok) 5788 end if 5789 end do 5790 end do 5791 end do 5792 end do 5793 end do 5794 end do 5795 5796 end subroutine ddb_get_d3matr
m_ddb/ddb_get_dchidet [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_dchidet
FUNCTION
Reads the non-linear optical susceptibility tensor and the first-order change in the linear dielectric susceptibility
INPUTS
ddb<type(ddb_type)>=Derivative Database. ramansr= if /= 0, impose sum rule on first-order derivatives of the electronic susceptibility with respect to atomic displacements nlflag= if =3, only the non-linear optical susceptibilities is computed
OUTPUT
dchide(3,3,3) = non-linear optical coefficients dchidt(natom,3,3,3) = first-order change of the electronic dielectric tensor induced by an individual atomic displacement iblok=Index of the block containing the data. 0 if block is not found. The caller should check the returned value.
SOURCE
3827 integer function ddb_get_dchidet(ddb, ramansr, nlflag, dchide, dchidt) result(iblok) 3828 3829 !Arguments ------------------------------- 3830 !scalars 3831 integer,intent(in) :: ramansr, nlflag 3832 class(ddb_type),intent(in) :: ddb 3833 !arrays 3834 real(dp),intent(out) :: dchide(3,3,3),dchidt(ddb%natom,3,3,3) 3835 3836 !Local variables ------------------------- 3837 !scalars 3838 integer :: rftyp 3839 !arrays 3840 integer :: rfelfd(4),rfphon(4),rfstrs(4) 3841 real(dp) :: qphnrm(3),qphon(3,3) 3842 3843 ! ********************************************************************* 3844 3845 qphon(:,:) = zero 3846 qphnrm(:) = one 3847 ! rfphon(1) = 1 ; rfphon(2:3) = 0 3848 rfelfd(:) = 2 3849 rfstrs(:) = 0 3850 rftyp = BLKTYP_d3E_xx 3851 3852 if (nlflag < 3) then 3853 rfphon(1) = 1 ; rfphon(2:3) = 0 3854 else 3855 rfphon(1) = 0 ; rfphon(2:3) = 0 3856 end if 3857 3858 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 3859 3860 if (iblok /= 0) then 3861 call dtchi(ddb%val(:,:,iblok),dchide,dchidt,ddb%mpert,ddb%natom,ramansr,nlflag) 3862 else 3863 ! Let the caller handle the error. 3864 dchide = huge(one); dchidt = huge(one) 3865 end if 3866 3867 end function ddb_get_dchidet
m_ddb/ddb_get_dielt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_dielt
FUNCTION
Reads the electronic dielectric tensor from the DDB file
INPUTS
ddb<type(ddb_type)>=Derivative database. rftyp = 1 if non-stationary block 2 if stationary block 3 if third order derivatives
OUTPUT
dielt(3,3) = Macroscopic dielectric tensor (electronic contribution) iblok=Index of the block containing the data. 0 if block is not found.
NOTES
dielt is initialized to one_3D if the derivatives are not available in the DDB file.
SOURCE
3660 integer function ddb_get_dielt(ddb, rftyp, dielt) result(iblok) 3661 3662 !Arguments ------------------------------- 3663 !scalars 3664 integer,intent(in) :: rftyp 3665 class(ddb_type),intent(in) :: ddb 3666 !arrays 3667 real(dp),intent(out) :: dielt(3,3) 3668 3669 !Local variables ------------------------- 3670 !scalars 3671 integer :: mpert 3672 character(len=1000) :: msg 3673 !arrays 3674 integer :: rfelfd(4),rfphon(4),rfstrs(4) 3675 real(dp) :: qphnrm(3),qphon(3,3) 3676 real(dp),allocatable :: tmpval(:,:,:,:) 3677 3678 ! ********************************************************************* 3679 3680 ! Look for the Gamma Block in the DDB 3681 qphon(:,1)=zero 3682 qphnrm(1)=zero 3683 rfphon(1:2)=0 3684 rfelfd(1:2)=2 3685 rfstrs(1:2)=0 3686 3687 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 3688 3689 ! Read the dielectric tensor only if the Gamma-block was found in the DDB 3690 ! In case it was not found, iblok = 0 3691 dielt=zero; dielt(1,1)=one; dielt(2,2)=one; dielt(3,3)=one 3692 3693 if (iblok/=0) then 3694 !Extration of dielectric tensor 3695 mpert = ddb%mpert 3696 3697 ABI_MALLOC(tmpval,(3,mpert,3,mpert)) 3698 tmpval(:,:,:,:) = reshape(ddb%val(1,:,iblok), shape = (/3,mpert,3,mpert/)) 3699 dielt=tmpval(1:3,ddb%natom+2,1:3,ddb%natom+2) 3700 3701 write(msg,'(a,3es16.6,3es16.6,3es16.6)' )& 3702 ' Dielectric Tensor ',& 3703 dielt(1,1),dielt(1,2),dielt(1,3),& 3704 dielt(2,1),dielt(2,2),dielt(2,3),& 3705 dielt(3,1),dielt(3,2),dielt(3,3) 3706 3707 call wrtout(std_out,msg) 3708 3709 ABI_FREE(tmpval) 3710 end if ! iblok not found 3711 3712 end function ddb_get_dielt
m_ddb/ddb_get_dielt_zeff [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_dielt_zeff
FUNCTION
Reads the Dielectric Tensor and the Effective Charges from the DDB file Impose the charge neutrality on the effective charges and eventually select some parts of the effective charges
INPUTS
Crystal<type(crystal_t)>=Crystal structure parameters rftyp = 1 if non-stationary block 2 if stationary block 3 if third order derivatives chneut=(0 => no ASR, 1 => equal repartition,2 => weighted repartition ) selectz=selection of some parts of the effective charge tensor attached to one atom. (0=> no selection, 1=> trace only, 2=> symmetric part only)
SIDE EFFECTS
ddb<type(ddb_type)>= The block with the effective charges is modified if charge neutrality is imposed.
OUTPUT
iblok=Index of the block containing the data. 0 if block is not found. dielt(3,3) = Macroscopic dielectric tensor zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement [zeff_raw(3,3,natom)]=effective charge on each atom before enforcing charge-neutrality.
NOTES
dielt and zeff are initialized to one_3D and zero if the derivatives are not available in the DDB file.
SOURCE
3572 integer function ddb_get_dielt_zeff(ddb, crystal, rftyp, chneut, selectz, dielt, zeff, zeff_raw) result(iblok) 3573 3574 !Arguments ------------------------------- 3575 !scalars 3576 integer,intent(in) :: rftyp,chneut,selectz 3577 class(ddb_type),intent(inout) :: ddb 3578 type(crystal_t),intent(in) :: crystal 3579 !arrays 3580 real(dp),intent(out) :: dielt(3,3),zeff(3,3,crystal%natom) 3581 real(dp),optional,intent(out) :: zeff_raw(3,3,crystal%natom) 3582 3583 !Local variables ------------------------- 3584 !scalars 3585 integer :: ii 3586 character(len=500) :: msg 3587 !arrays 3588 integer :: rfelfd(4),rfphon(4),rfstrs(4) 3589 real(dp) :: qphnrm(3),qphon(3,3), my_zeff_raw(3,3,crystal%natom) 3590 3591 ! ********************************************************************* 3592 3593 ! Look for the Gamma Block in the DDB 3594 qphon(:,1)=zero 3595 qphnrm(1)=zero 3596 rfphon(1:2)=1 3597 rfelfd(1:2)=2 3598 rfstrs(1:2)=0 3599 3600 !write(std_out,*)"ddb%mpert",ddb%mpert 3601 call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, rftyp) 3602 3603 ! Compute effective charges and dielectric tensor only if the Gamma-blok was found in the DDB 3604 ! In case it was not found, iblok = 0 3605 zeff=zero; dielt=zero; dielt(1,1)=one; dielt(2,2)=one; dielt(3,3)=one 3606 my_zeff_raw = zero 3607 3608 if (iblok /= 0) then 3609 write(msg, '(2a,(80a),4a)' ) ch10,('=',ii=1,80),ch10,ch10,& 3610 ' Dielectric Tensor and Effective Charges ',ch10 3611 call wrtout([std_out, ab_out], msg) 3612 3613 ! Make the imaginary part of the Gamma block vanish 3614 write(msg, '(5a)' ) ch10,& 3615 ' anaddb : Zero the imaginary part of the Dynamical Matrix at Gamma,',ch10,& 3616 ' and impose the ASR on the effective charges ',ch10 3617 call wrtout([std_out, ab_out], msg) 3618 3619 ! Extrac Zeff before enforcing sum rule. 3620 call dtech9(ddb%val, dielt, iblok, ddb%mpert, ddb%natom, ddb%nblok, my_zeff_raw, unit=dev_null) 3621 3622 ! Impose the charge neutrality on the effective charges and eventually select some parts of the effective charges 3623 call chneu9(chneut,ddb%val(:,:,iblok),ddb%mpert,ddb%natom,ddb%ntypat,selectz,Crystal%typat,Crystal%zion) 3624 3625 ! Extraction of the dielectric tensor and the effective charges 3626 call dtech9(ddb%val, dielt, iblok, ddb%mpert, ddb%natom, ddb%nblok, zeff) 3627 3628 end if ! iblok not found 3629 3630 if (present(zeff_raw)) zeff_raw = my_zeff_raw 3631 3632 end function ddb_get_dielt_zeff
m_ddb/ddb_get_etotal [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_etotal
FUNCTION
Read the GS total energy from the DDB file.
INPUTS
ddb<type(ddb_type)>=Derivative database.
OUTPUT
etotal=GS Total energy in Hartree iblok=Index of the block in the DDB file. 0 if not found.
SOURCE
3502 integer function ddb_get_etotal(ddb, etotal) result(iblok) 3503 3504 !Arguments ------------------------------- 3505 !scalars 3506 real(dp),intent(out) :: etotal 3507 class(ddb_type),intent(in) :: ddb 3508 3509 !Local variables ------------------------- 3510 !scalars 3511 integer :: rftyp 3512 !arrays 3513 integer :: rfelfd(4),rfphon(4),rfstrs(4) 3514 real(dp) :: qphnrm(3),qphon(3,3) 3515 3516 ! ********************************************************************* 3517 3518 ! Extract the block with the total energy 3519 qphon(:,:) = zero 3520 qphnrm(:) = zero 3521 rfphon(:) = 0 3522 rfelfd(:) = 0 3523 rfstrs(:) = 0 3524 rftyp = BLKTYP_d0E_xx 3525 3526 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 3527 3528 if (iblok /= 0) then 3529 etotal = ddb%val(1,1,iblok) 3530 else 3531 etotal = huge(one) 3532 end if 3533 3534 end function ddb_get_etotal
m_ddb/ddb_get_gred [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_gred
FUNCTION
Get the forces in reduced coordinates (Hartree).
INPUTS
ddb<type(ddb_type)>=Derivative database. relaxat 0 => without relaxation of the atoms 1 => with relaxation of the atoms relaxstr 0 => without relaxed lattice constants 1 => with relaxed lattice constants at constrained polarization
OUTPUT
gred(3,natom)=the gradient of the total energy with respect to change of reduced coordinates iblok=Index of the block containing the data. 0 if block is not found.
NOTES
SOURCE
3960 integer function ddb_get_gred(ddb, gred, relaxat, relaxstr) result(iblok) 3961 3962 !Arguments ------------------------------- 3963 !scalars 3964 class(ddb_type),intent(in) :: ddb 3965 integer, intent(in) :: relaxat 3966 integer, intent(in) :: relaxstr 3967 !arrays 3968 real(dp),intent(out),allocatable :: gred(:,:) 3969 3970 !Local variables ------------------------- 3971 !scalars 3972 integer :: natom 3973 integer :: idir, iatom, index 3974 !arrays 3975 integer :: rfelfd(4),rfphon(4),rfstrs(4) 3976 real(dp) :: qphnrm(3),qphon(3,3) 3977 3978 ! ********************************************************************* 3979 3980 natom = ddb%natom 3981 3982 ABI_MALLOC(gred, (3, natom)) 3983 3984 qphon(:,:) = zero; qphnrm(:) = zero 3985 rfphon(:) = 0; rfstrs(:) = 0; rfelfd(:) = 2 3986 if (relaxat == 1) rfphon(:) = 1 3987 if (relaxstr == 1) rfstrs(:) = 3 3988 3989 ! GA: I dont see why relaxstr is relevant as an input 3990 call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, BLKTYP_d1E_xx) 3991 3992 if (iblok/=0) then 3993 if (relaxat == 1) then 3994 index = 0 3995 do iatom = 1, natom 3996 do idir = 1, 3 3997 index = index+1 3998 gred(idir, iatom) = ddb%val(1, index, iblok) 3999 end do 4000 end do 4001 end if 4002 end if 4003 4004 end function ddb_get_gred
m_ddb/ddb_get_pel [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_pel
FUNCTION
Get the electronic polarizability vector from the database.
INPUTS
ddb<type(ddb_type)>=Derivative database. relaxat 0 => without relaxation of the atoms 1 => with relaxation of the atoms relaxstr 0 => without relaxed lattice constants 1 => with relaxed lattice constants at constrained polarization
OUTPUT
pel(3) = Macroscopic polarizability vector (electronic contribution) iblok=Index of the block containing the data. 0 if block is not found.
NOTES
SOURCE
3897 integer function ddb_get_pel(ddb, pel, relaxat, relaxstr) result(iblok) 3898 3899 !Arguments ------------------------------- 3900 !scalars 3901 class(ddb_type),intent(in) :: ddb 3902 integer, intent(in) :: relaxat 3903 integer, intent(in) :: relaxstr 3904 !arrays 3905 real(dp),intent(out) :: pel(3) 3906 3907 !Local variables ------------------------- 3908 !scalars 3909 integer :: natom 3910 !arrays 3911 integer :: rfelfd(4),rfphon(4),rfstrs(4) 3912 real(dp) :: qphnrm(3),qphon(3,3) 3913 3914 ! ********************************************************************* 3915 3916 natom = ddb%natom 3917 3918 qphon(:,:) = zero; qphnrm(:) = zero 3919 rfphon(:) = 0; rfstrs(:) = 0; rfelfd(:) = 2 3920 if (relaxat == 1) rfphon(:) = 1 3921 if (relaxstr == 1) rfstrs(:) = 3 3922 3923 call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, BLKTYP_d1E_xx) 3924 3925 if (iblok/=0) then 3926 pel(1:3) = ddb%val(1, 3*natom+4:3*natom+6, iblok) 3927 end if 3928 3929 end function ddb_get_pel
m_ddb/ddb_get_quadrupoles [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_quadrupoles
FUNCTION
Reads the Dynamic Quadrupoles or the P^(1) tensor from the DDB file
INPUTS
ddb<type(ddb_type)>=Derivative database. ddb_version = 6 digit integer giving date. To mantain compatibility with old DDB files. lwsym = 0 do not symmetrize the tensor wrt efield and qvec derivative |-> 1st gradient of polarization response to atomic displacement = 1 symmetrize the tensor wrt efield and qvec derivative |-> dynamic quadrupoles rftyp = 1 if non-stationary block 2 if stationary block 3 if third order derivatives 33 if long wave third order derivatives
OUTPUT
quadrupoles(3,3,3,natom) = Dynamic Quadrupole tensor iblok=Index of the block containing the data. 0 if block is not found.
NOTES
quadrupoles is initialized to zero if the derivatives are not available in the DDB file.
SOURCE
3746 integer function ddb_get_quadrupoles(ddb, ddb_version, lwsym, rftyp, quadrupoles) result(iblok) 3747 3748 !Arguments ------------------------------- 3749 !scalars 3750 integer,intent(in) :: ddb_version,lwsym,rftyp 3751 class(ddb_type),intent(in) :: ddb 3752 !arrays 3753 real(dp),intent(out) :: quadrupoles(3,3,3,ddb%natom) 3754 3755 !Local variables ------------------------- 3756 !scalars 3757 character(len=500) :: msg 3758 !arrays 3759 integer :: rfelfd(4),rfphon(4),rfstrs(4) 3760 integer :: rfqvec(4) 3761 real(dp) :: qphnrm(3),qphon(3,3) 3762 3763 ! ********************************************************************* 3764 3765 ! Look for the Gamma Block in the DDB 3766 qphon(:,:)=zero 3767 qphnrm(:)=one 3768 rfphon(:)=0 3769 rfelfd(:)=0 3770 rfqvec(:)=0 3771 rfstrs(:)=0 3772 rfelfd(1)=2 3773 rfphon(2)=1 3774 rfqvec(3)=1 3775 3776 call ddb%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,rftyp,rfqvec=rfqvec) 3777 3778 ! Compute the quadrupole tensor only if the Gamma-block was found in the DDB 3779 ! In case it was not found, iblok = 0 3780 quadrupoles=zero 3781 3782 if (iblok /= 0) then 3783 write(msg,'(2a)') ch10, ' Extract quadrupoles or P^(1) coefficients from 3DTE' 3784 call wrtout(std_out, msg) 3785 3786 if (lwsym==1) then 3787 write(msg, '(3a)' ) ch10, ' Dynamical Quadrupoles Tensor (units: e Bohr)',ch10 3788 else if (lwsym==0) then 3789 write(msg, '(3a)' ) ch10, & 3790 ' First moment of Polarization induced by atomic displacement (1/ucvol factor not included) (units: e Bohr) ',ch10 3791 endif 3792 call wrtout([std_out, ab_out], msg) 3793 3794 call dtqdrp(ddb%val(:,:,iblok),ddb_version,lwsym,ddb%mpert,ddb%natom,quadrupoles) 3795 end if 3796 3797 end function ddb_get_quadrupoles
m_ddb/ddb_get_strten [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_get_strten
FUNCTION
Get the stress tensor.
INPUTS
ddb<type(ddb_type)>=Derivative database. relaxat 0 => without relaxation of the atoms 1 => with relaxation of the atoms relaxstr 0 => without relaxed lattice constants 1 => with relaxed lattice constants at constrained polarization
OUTPUT
strten(6)=the stress tensor in cartesian coordinates. iblok=Index of the block containing the data. 0 if block is not found.
NOTES
SOURCE
4034 integer function ddb_get_strten(ddb, strten, relaxat, relaxstr) result(iblok) 4035 4036 !Arguments ------------------------------- 4037 !scalars 4038 class(ddb_type),intent(in) :: ddb 4039 integer, intent(in) :: relaxat 4040 integer, intent(in) :: relaxstr 4041 !arrays 4042 real(dp),intent(out) :: strten(6) 4043 4044 !Local variables ------------------------- 4045 !scalars 4046 integer :: natom 4047 integer :: ii, index 4048 !arrays 4049 integer :: rfelfd(4),rfphon(4),rfstrs(4) 4050 real(dp) :: qphnrm(3),qphon(3,3) 4051 4052 ! ********************************************************************* 4053 4054 natom = ddb%natom 4055 4056 qphon(:,:) = zero; qphnrm(:) = zero 4057 rfphon(:) = 0; rfstrs(:) = 0; rfelfd(:) = 2 4058 if (relaxat == 1) rfphon(:) = 1 4059 if (relaxstr == 1) rfstrs(:) = 3 4060 4061 ! GA: I dont see why relaxstr is relevant as an input 4062 call ddb%get_block(iblok, qphon, qphnrm, rfphon, rfelfd, rfstrs, BLKTYP_d1E_xx) 4063 4064 if (iblok/=0) then 4065 if (relaxstr == 1) then 4066 index = 3*natom+6 4067 do ii = 1, 6 4068 index = index+1 4069 strten(ii) = ddb%val(1, index, iblok) 4070 end do 4071 end if 4072 end if 4073 4074 end function ddb_get_strten
m_ddb/ddb_init [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_init
FUNCTION
Initialize a new ddb object for the current calculation.
INPUTS
ddb=the new ddb object
dtset=dtset object of the current calculation
nblok=number of blocks
mpert=maximum number of perturbations (atom displacements + electric field + ...)
with_d0E=this ddb contains 0th order derivatives
with_d1E=this ddb contains 1st order derivatives
with_d2E=this ddb contains 2nd order derivatives
with_d3E=this ddb contains 3rd order derivatives
with_d2eig=this ddb contains 2nd order derivatives of eigenvalues
mband='number of bands' dimension of the d2eig array.
Should actually correspond to the maximum number of bands for one kpoint
multiplied by the number of spin polarization (mband*nsppol).
nkpt=number of kpoints
kpt=reduced coordinates of kpoints
SOURCE
447 subroutine ddb_init(ddb, dtset, nblok, mpert, & 448 mband, nkpt, kpt,& 449 with_d0E, with_d1E, with_d2E, with_d3E, with_d2eig) 450 451 !Arguments ------------------------------- 452 class(ddb_type),intent(inout) :: ddb 453 type(dataset_type),intent(in) :: dtset 454 integer,intent(in) :: nblok, mpert 455 integer,intent(in),optional :: mband,nkpt 456 real(dp),intent(in),optional :: kpt(:,:) 457 logical,intent(in),optional :: with_d0E, with_d1E, with_d2E, with_d3E, with_d2eig 458 459 !Local variables ------------------------------- 460 integer :: msize_, ii, ikpt 461 logical :: with_d0E_, with_d1E_, with_d2E_, with_d3E_, with_d2eig_ 462 463 ! ************************************************************************ 464 465 with_d0E_ = .false. ; if (present(with_d0E)) with_d0E_ = with_d0E 466 with_d1E_ = .false. ; if (present(with_d1E)) with_d1E_ = with_d1E 467 with_d2E_ = .false. ; if (present(with_d2E)) with_d2E_ = with_d2E 468 with_d3E_ = .false. ; if (present(with_d3E)) with_d3E_ = with_d3E 469 with_d2eig_ = .false. ; if (present(with_d2eig)) with_d2eig_ = with_d2eig 470 471 msize_ = 0 472 if (with_d0E_) msize_ = 1 473 if (with_d1E_) msize_ = 3 * mpert 474 if (with_d2E_ .or. with_d2eig_) msize_ = 3 * mpert * 3 * mpert 475 if (with_d3E_) msize_ = 3 * mpert * 3 * mpert * 3 * mpert 476 477 call ddb%malloc(msize_, nblok, dtset%natom, dtset%ntypat, mpert) 478 479 ddb%occopt = dtset%occopt 480 ddb%prtvol = dtset%prtvol 481 482 ddb%rprim(:,:) = dtset%rprim_orig(1:3,1:3,1) 483 ddb%acell(:) = dtset%acell_orig(1:3,1) 484 485 call matr3inv(ddb%rprim, ddb%gprim) 486 487 ddb%qpt(:,:) = zero 488 ddb%nrm(:,:) = one 489 if (with_d0E_) then 490 ddb%typ(:) = BLKTYP_d0E_xx 491 else if (with_d1E_) then 492 ddb%typ(:) = BLKTYP_d1E_xx 493 else if (with_d2E_) then 494 ddb%typ(:) = BLKTYP_d2E_ns 495 else if (with_d3E_) then 496 ddb%typ(:) = BLKTYP_d3E_xx 497 else if (with_d2eig_) then 498 ddb%typ(:) = BLKTYP_d2eig_re 499 end if 500 501 ddb%flg(:,:) = 0 502 ddb%amu(:) = dtset%amu_orig(:,1) 503 504 ddb%nsppol = dtset%nsppol 505 506 if (present(mband)) then 507 ddb%nband = mband 508 else 509 ddb%nband = dtset%mband * ddb%nsppol 510 end if 511 512 if (present(nkpt)) then 513 ddb%nkpt = nkpt 514 else 515 ddb%nkpt = dtset%nkpt 516 end if 517 518 ! TODO: Allocate d2eig here instead of leaving it to the calling routine. 519 if (with_d2eig_) then 520 call ddb%malloc_d2eig(ddb%nband, ddb%nkpt) 521 end if 522 523 if (present(kpt)) then 524 do ikpt=1,ddb%nkpt 525 do ii = 1,3 526 ddb%kpt(ii,ikpt) = kpt(ii,ikpt) 527 end do 528 end do 529 end if 530 531 end subroutine ddb_init
m_ddb/ddb_lw_copy [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_lw_copy
FUNCTION
Copy the ddb object after reading the long wave 3rd order derivatives into a new ddb_lw and resizes ddb as for 2nd order derivatives
INPUTS
ddb (INOUT) = ddb block datastructure mpert =maximum number of ipert natom= number of atoms in unit cell ntypat= number of atom types
OUTPUT
ddb_lw= ddb block datastructure
SOURCE
6699 subroutine ddb_lw_copy(ddb,ddb_lw,mpert,natom,ntypat) 6700 6701 !Arguments ------------------------------- 6702 !scalars 6703 integer,intent(in) :: mpert,natom,ntypat 6704 !arrays 6705 type(ddb_type),intent(inout) :: ddb 6706 type(ddb_type),intent(out) :: ddb_lw 6707 6708 !Local variables ------------------------- 6709 !scalars 6710 integer :: ii,nblok,nsize,cnt 6711 6712 ! ********************************************************************* 6713 6714 call ddb%copy(ddb_lw) 6715 call ddb%free() 6716 nsize=3*mpert*3*mpert 6717 nblok=ddb_lw%nblok-count(ddb_lw%typ(:)==BLKTYP_d3E_lw) 6718 call ddb%malloc(nsize, nblok, natom, ntypat, mpert) 6719 6720 ! Copy dimensions and static variables. 6721 ddb%msize = nsize 6722 ddb%mpert = ddb_lw%mpert 6723 ddb%nblok = nblok 6724 ddb%natom = ddb_lw%natom 6725 ddb%ntypat = ddb_lw%ntypat 6726 ddb%occopt = ddb_lw%occopt 6727 ddb%prtvol = ddb_lw%prtvol 6728 6729 ddb%rprim = ddb_lw%rprim 6730 ddb%gprim = ddb_lw%gprim 6731 ddb%acell = ddb_lw%acell 6732 6733 ! Copy the allocatable arrays. 6734 ddb%amu(:) = ddb_lw%amu(:) 6735 cnt = 0 6736 do ii=1,ddb_lw%nblok 6737 if (ddb_lw%typ(ii)/=BLKTYP_d3E_lw) then 6738 cnt = cnt + 1 6739 ddb%flg(:,cnt) = ddb_lw%flg(1:nsize,ii) 6740 ddb%val(:,:,cnt) = ddb_lw%val(:,1:nsize,ii) 6741 ddb%typ(cnt) = ddb_lw%typ(ii) 6742 ddb%nrm(:,cnt) = ddb_lw%nrm(:,ii) 6743 ddb%qpt(:,cnt) = ddb_lw%qpt(:,ii) 6744 end if 6745 end do 6746 6747 end subroutine ddb_lw_copy
m_ddb/ddb_malloc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_malloc
FUNCTION
Allocate dynamic memory.
INPUTS
msize=maximum size of one block of the ddb
(e.g. 3*mpert * 3*mpert)
nblok=number of blocks in the ddb
natom=number of atoms
ntypat=number of atom types
mpert=maximum number of perturbations
(atom displacements + electric field + ...)
nkpt=number of k-points. Optional, indicates the use of eig2d.
nband='number of bands' dimension of the d2eig array.
Should actually correspond to the maximum number of bands for one kpoint
multiplied by the number of spin polarization (mband*nsppol).
OUTPUT
SOURCE
637 subroutine ddb_malloc(ddb, msize, nblok, natom, ntypat, mpert, nkpt, nband) 638 639 !Arguments ------------------------------- 640 !array 641 class(ddb_type),intent(inout) :: ddb 642 integer,intent(in) :: msize,nblok,natom,ntypat,mpert 643 integer,intent(in),optional :: nkpt,nband 644 645 ! ************************************************************************ 646 647 ddb%msize = msize 648 ddb%nblok = nblok 649 ddb%natom = natom 650 !ddb%mpert = natom + MPERT_MAX 651 ddb%mpert = mpert 652 ddb%ntypat = ntypat 653 654 ! integer 655 ABI_CALLOC(ddb%flg, (msize, nblok)) 656 ABI_CALLOC(ddb%typ, (nblok)) 657 658 ! real 659 ABI_MALLOC(ddb%amu, (ntypat)) 660 ABI_MALLOC(ddb%nrm, (3, nblok)) 661 ABI_MALLOC(ddb%qpt, (9, nblok)) 662 ABI_MALLOC(ddb%val, (2, msize, nblok)) 663 ddb%val = huge(one) 664 665 ! FIXME: should really add nsppol argument (see thmeig). 666 if (present(nkpt) .and. present(nband)) then 667 call ddb%malloc_d2eig(nband, nkpt) 668 end if 669 670 end subroutine ddb_malloc
m_ddb/ddb_malloc_d2eig [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_malloc_d2eig
FUNCTION
Allocate dynamic memory for second derivatives of eigenvalues.
INPUTS
mband='number of bands' dimension of the d2eig array.
Should actually correspond to the maximum number of bands for one kpoint
multiplied by the number of spin polarization (mband*nsppol).
nkpt=number of kpoints
OUTPUT
SOURCE
692 subroutine ddb_malloc_d2eig(ddb, mband, nkpt) 693 694 !Arguments ------------------------------- 695 !array 696 class(ddb_type),intent(inout) :: ddb 697 integer,intent(in) :: mband, nkpt 698 699 ! ************************************************************************ 700 701 ddb%nband = mband / ddb%nsppol 702 ddb%nkpt = nkpt 703 ABI_MALLOC(ddb%kpt, (3, nkpt)) 704 ABI_MALLOC(ddb%eig2dval, (2, ddb%msize, mband, nkpt)) 705 706 end subroutine ddb_malloc_d2eig
m_ddb/ddb_merge_blocks [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_merge_blocks
FUNCTION
Merge block number iblok2 from ddb2 into block number iblok1 in ddb1.
INPUTS
ddb1=ddb object 1 ddb2=ddb object 2 iblok1=block index from ddb1 iblok2=block index from ddb2
OUTPUT
SOURCE
2868 subroutine ddb_merge_blocks(ddb1, ddb2, iblok1, iblok2) 2869 2870 !Arguments ------------------------------- 2871 !array 2872 class(ddb_type),intent(inout) :: ddb1 2873 type(ddb_type),intent(inout) :: ddb2 2874 integer,intent(in) :: iblok1 2875 integer,intent(in) :: iblok2 2876 2877 !local variables 2878 !scalars 2879 integer :: ii, blktyp, mpert1, mpert2 2880 integer :: idir1, idir2, idir3, ipert1, ipert2, ipert3 2881 real(dp),parameter :: qtol=2.0d-8 2882 !arrays 2883 real(dp), allocatable :: d1matr(:,:,:) 2884 real(dp), allocatable :: d2matr(:,:,:,:,:) 2885 real(dp), allocatable :: d3matr(:,:,:,:,:,:,:) 2886 integer, allocatable :: d1flg(:,:) 2887 integer, allocatable :: d2flg(:,:,:,:) 2888 integer, allocatable :: d3flg(:,:,:,:,:,:) 2889 2890 ! ************************************************************************ 2891 2892 ! Note that ddb and ddb2 may have a different values for mpert 2893 !mpert = min(ddb1%mpert, ddb2%mpert) 2894 mpert1 = ddb1%mpert 2895 mpert2 = ddb2%mpert 2896 2897 ! Add the blok to the output ddb 2898 blktyp = ddb2%typ(iblok2) 2899 ddb1%typ(iblok1) = blktyp 2900 2901 ! Copy q-point 2902 do ii=1,9 2903 ddb1%qpt(ii,iblok1) = ddb2%qpt(ii,iblok2) 2904 end do 2905 do ii=1,3 2906 ddb1%nrm(ii,iblok1) = ddb2%nrm(ii,iblok2) 2907 end do 2908 2909 if (is_type_d0E(blktyp)) then 2910 ! -------------- 2911 ! Copy d0E block 2912 ! -------------- 2913 if (ddb2%flg(1,iblok2) > 0) then 2914 ddb1%val(1,1,iblok1) = ddb2%val(1,1,iblok2) 2915 ddb1%val(2,1,iblok1) = ddb2%val(2,1,iblok2) 2916 ddb1%flg(1,iblok1) = ddb2%flg(1,iblok2) 2917 end if 2918 2919 else if (is_type_d1E(blktyp)) then 2920 ! -------------- 2921 ! Copy d1E block 2922 ! -------------- 2923 call ddb2%get_d1matr(iblok2, d1matr, d1flg) 2924 !call ddb%set_d1matr(iblok, d1matr, d1flg) 2925 ii=0 2926 do ipert1=1,mpert1 2927 do idir1=1,3 2928 ii=ii+1 2929 if (ipert1 <= mpert2) then 2930 if (d1flg(idir1,ipert1)>0) then 2931 ddb1%val(1,ii,iblok1) = d1matr(1,idir1,ipert1) 2932 ddb1%val(2,ii,iblok1) = d1matr(2,idir1,ipert1) 2933 ddb1%flg(ii,iblok1) = d1flg(idir1,ipert1) 2934 end if 2935 end if 2936 end do 2937 end do 2938 ABI_SFREE(d1matr) 2939 ABI_SFREE(d1flg) 2940 2941 else if (is_type_d2E(blktyp)) then 2942 ! -------------- 2943 ! Copy d2E block 2944 ! -------------- 2945 call ddb2%get_d2matr(iblok2, d2matr, d2flg) 2946 !call ddb%set_d2matr(iblok, d2matr, d2flg) 2947 ii=0 2948 do ipert2=1,mpert1 2949 do idir2=1,3 2950 do ipert1=1,mpert1 2951 do idir1=1,3 2952 ii=ii+1 2953 if ((ipert1 <= mpert2).and.(ipert2<=mpert2)) then 2954 if (d2flg(idir1,ipert1,idir2,ipert2)>0) then 2955 ddb1%val(1,ii,iblok1) = d2matr(1,idir1,ipert1,idir2,ipert2) 2956 ddb1%val(2,ii,iblok1) = d2matr(2,idir1,ipert1,idir2,ipert2) 2957 ddb1%flg(ii,iblok1) = d2flg(idir1,ipert1,idir2,ipert2) 2958 end if 2959 end if 2960 end do 2961 end do 2962 end do 2963 end do 2964 ABI_SFREE(d2matr) 2965 ABI_SFREE(d2flg) 2966 2967 else if (is_type_d3E(blktyp)) then 2968 ! -------------- 2969 ! Copy d3E block 2970 ! -------------- 2971 call ddb2%get_d3matr(iblok2, d3matr, d3flg) 2972 !call ddb%set_d3matr(iblok, d3matr, d3flg) 2973 ii=0 2974 do ipert1=1,mpert1 2975 do idir1=1,3 2976 do ipert2=1,mpert1 2977 do idir2=1,3 2978 do ipert3=1,mpert1 2979 do idir3=1,3 2980 2981 2982 !ii=ii+1 ! GA: This is not equivalent 2983 ii = idir1 + 3*((ipert1-1)+mpert1*((idir2-1) & 2984 + 3*((ipert2-1)+mpert1*((idir3-1) & 2985 + 3*(ipert3-1))))) 2986 2987 ! Note that the loop order (1,2,3) is not really consistent 2988 ! with the d2E case (2, 1) 2989 ! TODO Clean this up 2990 2991 if ((ipert1 <= mpert2).and.(ipert2<=mpert2).and.(ipert3<=mpert2)) then 2992 if (d3flg(idir1,ipert1,idir2,ipert2,idir3,ipert3)>0) then 2993 ddb1%flg(ii,iblok1) = d3flg(idir1,ipert1,idir2,ipert2,idir3,ipert3) 2994 ddb1%val(:,ii,iblok1) = d3matr(:,idir1,ipert1,idir2,ipert2,idir3,ipert3) 2995 end if 2996 end if 2997 2998 end do 2999 end do 3000 end do 3001 end do 3002 end do 3003 end do 3004 ABI_SFREE(d3matr) 3005 ABI_SFREE(d3flg) 3006 3007 else if (is_type_d2eig(blktyp)) then 3008 ! ---------------- 3009 ! Skip d2eig block 3010 ! ---------------- 3011 3012 ! TODO need a function ddb_merge_d2eig(filename1, filename2, ) 3013 3014 end if ! blktyp 3015 3016 end subroutine ddb_merge_blocks
m_ddb/ddb_read_block_txt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_block_txt
FUNCTION
Read the next block of data from a DDB in text format.
INPUTS
iblok=the blok index to be assigned mpert=maximum number of ipert msize=maximum size of the arrays flags and values nunit=unit number for the data block file
OUTPUT
(see side effects)
SIDE EFFECTS
Input/Output ddb = ddb block datastructure ddb%typ=type of the block: 0 => total energy 1 => second-order energy derivatives, non-stationary block 2 => second-order energy derivatives, stationary block 3 => third-order energy derivatives 4 => first-order energy derivatives: forces, stresses and polarization 5 => second-order eigenvalue derivatives ddb%flg(msize)=flag for every matrix element (0=> the element is not in the data block), (1=> the element is in the data blok) ddb%qpt(9)=wavevector of the perturbation(s). The elements from 1 to 3 are used if we are dealing with the 2nd derivative of total energy (only one wavevector), while all elements are used in case of a third order derivative of total energy (three wavevector could be present) ddb%nrm(3)=normalization factors for the three allowed wavevectors. ddb%val(2,msize)=real(dp), complex, value of the matrix elements that are present in the data block [blkval2(2,msize,mband,nkpt)]= value of the matrix elements that are present in a block of EIGR2D/EIGI2D
NOTES
only executed by one processor.
SOURCE
1659 subroutine ddb_read_block_txt(ddb,iblok,mband,mpert,msize,nkpt,nunit,& 1660 blkval2,kpt) !optional 1661 1662 !Arguments ------------------------------- 1663 !scalars 1664 integer,intent(in) :: mband,mpert,msize,nkpt,nunit 1665 integer, intent(in) :: iblok 1666 !logical, intent(in), optional :: eig2d 1667 class(ddb_type),intent(inout) :: ddb 1668 !arrays 1669 real(dp),intent(out),optional :: kpt(3,nkpt) 1670 real(dp),intent(out),optional :: blkval2(2,msize,mband,nkpt) 1671 1672 !Local variables ------------------------- 1673 !scalars 1674 integer :: band,iband,idir1,idir2,idir3,ii,ikpt,index,ipert1,ipert2,ipert3,nelmts 1675 logical :: eig2d_ 1676 real(dp) :: ai,ar 1677 character(len=32) :: name 1678 character(len=500) :: msg 1679 1680 ! ********************************************************************* 1681 1682 ! Zero every flag 1683 ddb%flg(1:msize, iblok)=0 1684 1685 eig2d_ = .false. 1686 1687 1688 if(present(kpt).and.present(blkval2)) then 1689 ! GA: Weird that it is not allocated here 1690 blkval2(:,:,:,:)=zero 1691 kpt(:,:)=zero 1692 eig2d_ = .true. 1693 end if 1694 1695 !if (present(eig2d)) then 1696 ! eig2d_ = eig2d 1697 !end if 1698 1699 ! Read the block type and number of elements 1700 read(nunit,*) 1701 read(nunit, '(a32,12x,i12)' )name,nelmts 1702 1703 ! TODO: Replace with STRING_d2E, etc. 1704 ! GA: Note that older versions used the expression '2rd' instead of '2nd' 1705 ! So this substitution needs to be checked for backward compatibility. 1706 ! Also, the strings for d2eig and d2eig_brd are undistinguishable 1707 ! I don't think the d2eig_brd was ever read by abinit or anaddb. 1708 if(name==' 2nd derivatives (non-stat.) - ' .or. name==' 2rd derivatives (non-stat.) - ')then 1709 ddb%typ(iblok)=BLKTYP_d2E_ns 1710 else if(name==' 2nd derivatives (stationary) - ' .or. name==' 2rd derivatives (stationary) - ')then 1711 ddb%typ(iblok)=BLKTYP_d2E_st 1712 else if(name==' 3rd derivatives - ')then 1713 ddb%typ(iblok)=BLKTYP_d3E_xx 1714 else if(name==' Total energy - ')then 1715 ddb%typ(iblok)=BLKTYP_d0E_xx 1716 else if(name==' 1st derivatives - ')then 1717 ddb%typ(iblok)=BLKTYP_d1E_xx 1718 else if(name==' 2nd eigenvalue derivatives - ' .or. name==' 2rd eigenvalue derivatives - ')then 1719 ddb%typ(iblok)=BLKTYP_d2eig_re 1720 else if(name==' 3rd derivatives (long wave) - ')then 1721 ddb%typ(iblok)=BLKTYP_d3E_lw 1722 else if(name==' 2nd derivatives (MBC) - ')then 1723 ddb%typ(iblok)=BLKTYP_d2E_mbc 1724 else 1725 write(msg,'(6a)')& 1726 'The following string appears in the DDB in place of',& 1727 ' the block type description :',ch10,trim(name),ch10,& 1728 'Action: check your DDB.' 1729 ABI_ERROR(msg) 1730 end if 1731 1732 ! Read the 2nd derivative block 1733 if (is_type_d2E(ddb%typ(iblok))) then 1734 1735 ! First check if there is enough space to read it 1736 if(msize<(3*mpert*3*mpert))then 1737 write(msg,'(3a)')& 1738 'There is not enough space to read a second-derivative block.',ch10,& 1739 'Action: increase msize and recompile.' 1740 ABI_ERROR(msg) 1741 end if 1742 1743 ! Read the phonon wavevector 1744 read(nunit, '(4x,3es16.8,f6.1)' )(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok) 1745 1746 ! Read every element 1747 do ii=1,nelmts 1748 read(nunit,*)idir1,ipert1,idir2,ipert2,ar,ai 1749 index=idir1+3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1))) 1750 ddb%flg(index,iblok)=1 1751 ddb%val(1,index,iblok)=ar 1752 ddb%val(2,index,iblok)=ai 1753 end do 1754 1755 else if (is_type_d3E(ddb%typ(iblok))) then 1756 ! Read the 3rd derivative block 1757 1758 ! First check if there is enough space to read it 1759 if(msize<(3*mpert*3*mpert*3*mpert))then 1760 write(msg, '(a,a,a,i10,a,i10,a,a,a)' )& 1761 'There is not enough space to read a third-derivative block.',ch10,& 1762 'The size provided is only ',msize,' although ',3*mpert*3*mpert*3*mpert,' is needed.',ch10,& 1763 'Action: increase msize and recompile.' 1764 ABI_ERROR(msg) 1765 end if 1766 1767 ! Read the perturbation wavevectors 1768 read(nunit,'(4x,3es16.8,f6.1)')(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok) 1769 read(nunit,'(4x,3es16.8,f6.1)')(ddb%qpt(ii,iblok),ii=4,6),ddb%nrm(2,iblok) 1770 read(nunit,'(4x,3es16.8,f6.1)')(ddb%qpt(ii,iblok),ii=7,9),ddb%nrm(3,iblok) 1771 1772 ! Read every element 1773 do ii=1,nelmts 1774 read(nunit,*)idir1,ipert1,idir2,ipert2,idir3,ipert3,ar,ai 1775 index=idir1+ & 1776 3*((ipert1-1)+mpert*((idir2-1)+ & 1777 3*((ipert2-1)+mpert*((idir3-1)+3*(ipert3-1))))) 1778 ddb%flg(index,iblok)=1 1779 ddb%val(1,index,iblok)=ar 1780 ddb%val(2,index,iblok)=ai 1781 end do 1782 1783 1784 else if (is_type_d0E(ddb%typ(iblok))) then 1785 ! Read the total energy 1786 ! First check if there is enough space to read it 1787 if(msize<1)then 1788 write(msg, '(3a,i0,3a)' )& 1789 'There is not enough space to read a total energy block.',ch10,& 1790 'The size provided is only ',msize,' although 1 is needed.',ch10,& 1791 'Action: increase msize and recompile.' 1792 ABI_ERROR(msg) 1793 end if 1794 1795 ! Read the total energy 1796 read(nunit,'(2d22.14)')ar,ai 1797 ddb%flg(1,iblok)=1 1798 ddb%val(1,1,iblok)=ar 1799 ddb%val(2,1,iblok)=ai 1800 1801 1802 else if (is_type_d1E(ddb%typ(iblok))) then 1803 ! Read the 1st derivative block 1804 ! First check if there is enough space to read it 1805 if (msize < (3*mpert)) then 1806 write(msg, '(3a,i0,a,i0,3a)' )& 1807 'There is not enough space to read a first-derivative block.',ch10,& 1808 'The size provided is only ',msize,' although ',3*mpert,' is needed.',ch10,& 1809 'Action: increase msize and recompile.' 1810 ABI_ERROR(msg) 1811 end if 1812 1813 ! Read every element 1814 do ii=1,nelmts 1815 read(nunit,*)idir1,ipert1,ar,ai 1816 index=idir1 + 3*(ipert1 - 1) 1817 ddb%flg(index,iblok)=1 1818 ddb%val(1,index,iblok)=ar 1819 ddb%val(2,index,iblok)=ai 1820 end do 1821 1822 1823 else if (is_type_d2eig(ddb%typ(iblok))) then 1824 1825 ! Read the 2nd eigenvalue derivative block 1826 ! First check if there is enough space to read it 1827 if(msize<(3*mpert*3*mpert))then 1828 write(msg, '(3a,i0,a,i0,3a)' )& 1829 'There is not enough space to read a second-derivative block.',ch10,& 1830 'The size provided is only ',msize,' although ',3*mpert*3*mpert*mband*nkpt,' is needed.',ch10,& 1831 'Action: increase msize and recompile.' 1832 ABI_ERROR(msg) 1833 end if 1834 1835 ! Read the phonon wavevector 1836 read(nunit, '(4x,3es16.8,f6.1)' )(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok) 1837 1838 ! Read the K point and band 1839 if (eig2d_) then 1840 do ikpt=1,nkpt 1841 read(nunit, '(9x,3es16.8)')(kpt(ii,ikpt),ii=1,3) 1842 do iband=1,mband 1843 read(nunit, '(6x,i3)') band 1844 ! Read every element 1845 do ii=1,nelmts 1846 read(nunit,*)idir1,ipert1,idir2,ipert2,ar,ai 1847 index=idir1+3*((ipert1-1)+mpert*((idir2-1)+3*(ipert2-1))) 1848 ddb%flg(index,iblok)=1 1849 blkval2(1,index,iband,ikpt)=ar 1850 blkval2(2,index,iband,ikpt)=ai 1851 end do !nelmts 1852 end do !band 1853 end do !kpt 1854 end if 1855 end if 1856 1857 end subroutine ddb_read_block_txt
m_ddb/ddb_read_d0E_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_d0E_nc
FUNCTION
Read a DDB block containing 0th order derivatives of energy.
INPUTS
ncid=netcdf identifier of a file open in reading mode. iblok=index of the block we are setting. iblok_d0E=index of the block we are reading in the d0E group.
OUTPUT
SOURCE
5289 subroutine ddb_read_d0E_nc(ddb, ncid, iblok, iblok_d0E) 5290 5291 !Arguments ------------------------------- 5292 !scalars 5293 class(ddb_type),intent(inout) :: ddb 5294 integer,intent(in) :: ncid,iblok,iblok_d0E 5295 5296 !Local variables ------------------------- 5297 !scalars 5298 integer :: ncid_d0E 5299 integer :: ncerr 5300 !arrays 5301 integer :: flg(1) 5302 real(dp) :: val(1) 5303 5304 ! ************************************************************************ 5305 5306 ncid_d0E = nctk_idgroup(ncid, 'd0E') 5307 5308 ! Allocate temporary arrays 5309 ncerr = nf90_get_var(ncid_d0E, nctk_idname(ncid_d0E, 'matrix_values'), val, start=[1,1,iblok_d0E], count=[1,1,1]) 5310 NCF_CHECK(ncerr) 5311 ddb%val(1,1,iblok) = val(1) 5312 ddb%val(2,1,iblok) = zero 5313 ncerr = nf90_get_var(ncid_d0E, nctk_idname(ncid_d0E, 'matrix_mask'), flg, start=[1,iblok_d0E], count=[1,1]) 5314 NCF_CHECK(ncerr) 5315 ddb%flg(1,iblok) = flg(1) 5316 5317 end subroutine ddb_read_d0E_nc
m_ddb/ddb_read_d1E_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_d1E_nc
FUNCTION
Read a DDB block containing 1st order derivatives of energy.
INPUTS
ncid=netcdf identifier of a file open in reading mode. iblok=index of the block we are setting. iblok_d1E=index of the block we are reading in the d1E group.
OUTPUT
SOURCE
5339 subroutine ddb_read_d1E_nc(ddb, ncid, iblok, iblok_d1E) 5340 5341 !Arguments ------------------------------- 5342 !scalars 5343 class(ddb_type),intent(inout) :: ddb 5344 integer,intent(in) :: ncid,iblok,iblok_d1E 5345 5346 !Local variables ------------------------- 5347 !scalars 5348 integer :: ncid_d1E 5349 integer :: ncerr 5350 !arrays 5351 integer,allocatable :: flg_d1E(:,:) 5352 real(dp),allocatable :: matrix_d1E(:,:,:) 5353 5354 ! ************************************************************************ 5355 5356 ncid_d1E = nctk_idgroup(ncid, 'd1E') 5357 5358 ! Allocate temporary arrays 5359 ABI_MALLOC(matrix_d1E, (2,3,ddb%mpert)) 5360 ABI_MALLOC(flg_d1E, (3,ddb%mpert)) 5361 5362 ncerr = nf90_get_var(ncid_d1E, nctk_idname(ncid_d1E, 'matrix_values'), matrix_d1E, start=[1,1,1,iblok_d1E]) 5363 NCF_CHECK(ncerr) 5364 ncerr = nf90_get_var(ncid_d1E, nctk_idname(ncid_d1E, 'matrix_mask'), flg_d1E, start=[1,1,iblok_d1E]) 5365 NCF_CHECK(ncerr) 5366 5367 ! Reshape 5368 call ddb%set_d1matr(iblok, matrix_d1E, flg_d1E) 5369 5370 ! Free memory 5371 ABI_FREE(matrix_d1E) 5372 ABI_FREE(flg_d1E) 5373 5374 end subroutine ddb_read_d1E_nc
m_ddb/ddb_read_d2E_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_d2E_nc
FUNCTION
Read a DDB block containing 2nd order derivatives of energy.
INPUTS
ncid=netcdf identifier of a file open in reading mode. iblok=index of the block we are setting. iblok_d2E=index of the block we are reading in the d2E group.
OUTPUT
SOURCE
5395 subroutine ddb_read_d2E_nc(ddb, ncid, iblok, iblok_d2E) 5396 5397 !Arguments ------------------------------- 5398 !scalars 5399 class(ddb_type),intent(inout) :: ddb 5400 integer,intent(in) :: ncid,iblok,iblok_d2E 5401 5402 !Local variables ------------------------- 5403 !scalars 5404 integer :: ncid_d2E 5405 integer :: ncerr 5406 !arrays 5407 real(dp) :: qpt(3) 5408 integer,allocatable :: flg_d2E(:,:,:,:) 5409 real(dp),allocatable :: matrix_d2E(:,:,:,:,:) 5410 5411 ! ************************************************************************ 5412 5413 ncid_d2E = nctk_idgroup(ncid, 'd2E') 5414 5415 ! Allocate temporary arrays 5416 ABI_MALLOC(matrix_d2E, (2,3,ddb%mpert,3,ddb%mpert)) 5417 ABI_MALLOC(flg_d2E, (3,ddb%mpert,3,ddb%mpert)) 5418 5419 ncerr = nf90_get_var(ncid_d2E, nctk_idname(ncid_d2E, 'reduced_coordinates_of_qpoints'), qpt, start=[1,iblok_d2E]) 5420 NCF_CHECK(ncerr) 5421 ddb%qpt(1:3,iblok) = qpt(:) 5422 ncerr = nf90_get_var(ncid_d2E, nctk_idname(ncid_d2E, 'qpoints_normalization'), ddb%nrm(1,iblok), start=[iblok_d2E]) 5423 NCF_CHECK(ncerr) 5424 5425 ncerr = nf90_get_var(ncid_d2E, nctk_idname(ncid_d2E, 'matrix_values'), matrix_d2E, start=[1,1,1,1,1,iblok_d2E]) 5426 NCF_CHECK(ncerr) 5427 ncerr = nf90_get_var(ncid_d2E, nctk_idname(ncid_d2E, 'matrix_mask'), flg_d2E, start=[1,1,1,1,iblok_d2E]) 5428 NCF_CHECK(ncerr) 5429 5430 ! Reshape 5431 call ddb%set_d2matr(iblok, matrix_d2E, flg_d2E) 5432 5433 ! Free memory 5434 ABI_FREE(matrix_d2E) 5435 ABI_FREE(flg_d2E) 5436 5437 end subroutine ddb_read_d2E_nc
m_ddb/ddb_read_d2eig [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_d2eig
FUNCTION
Read the next DDB block containing second-order derivatives of eigenvalues and store it in block number iblok. The values of nband and nkpt must be set.
INPUTS
ddb_hdr=ddb header object with open file. iblok_store=the block index in the ddb object iblok_read=the block index in the ddb file
OUTPUT
NOTE The ddb object must be allocated becore calling this routine.
SOURCE
1885 subroutine ddb_read_d2eig(ddb, ddb_hdr, iblok_store, iblok_read, comm) 1886 1887 !Arguments ------------------------------- 1888 !scalars 1889 class(ddb_type),intent(inout) :: ddb 1890 type(ddb_hdr_type),intent(in) :: ddb_hdr 1891 integer, intent(in) :: iblok_store 1892 integer, intent(in),optional :: iblok_read 1893 integer, intent(in),optional :: comm 1894 1895 !Local variables ------------------------- 1896 !scalars 1897 integer,parameter :: master=0 1898 integer :: comm_ 1899 character(len=500) :: msg 1900 1901 ! ********************************************************************* 1902 1903 if (present(comm)) then 1904 comm_ = comm 1905 else 1906 comm_ = xmpi_comm_self 1907 end if 1908 1909 if (xmpi_comm_rank(comm_) == master) then 1910 1911 if (ddb_hdr%has_open_file_nc) then 1912 1913 ! Read the specified block and store it 1914 call ddb%read_d2eig_nc(ddb_hdr%ncid, iblok_store, iblok_read) 1915 1916 else if (ddb_hdr%has_open_file_txt) then 1917 1918 ! Read the next block and store it 1919 call ddb%read_d2eig_txt(ddb_hdr%unddb, iblok_store) 1920 1921 else 1922 write(msg, '(3a)' )& 1923 ! File has not beed open by ddb_hdr 1924 'Attempting to read from unopen file DDB.',ch10,& 1925 'Action: contact Abinit group.' 1926 ABI_ERROR(msg) 1927 end if 1928 1929 end if 1930 1931 ! GA: Should in principle broadcast the d2eig array, 1932 ! but I dont think it is required yet. 1933 1934 end subroutine ddb_read_d2eig
m_ddb/ddb_read_d2eig_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_d2eig_nc
FUNCTION
Read a DDB block containing 2nd order derivatives of eigenvalues.
INPUTS
ncid=netcdf identifier of a file open in reading mode. iblok=index of the block we are setting. iblok_d2eig=index of the block we are reading in the d2eig group.
OUTPUT
SOURCE
5534 subroutine ddb_read_d2eig_nc(ddb, ncid, iblok, iblok_d2eig) 5535 5536 !Arguments ------------------------------- 5537 !scalars 5538 class(ddb_type),intent(inout) :: ddb 5539 integer,intent(in) :: ncid,iblok 5540 integer,intent(in),optional :: iblok_d2eig 5541 5542 !Local variables ------------------------- 5543 !scalars 5544 integer :: ncid_d2eig,ncerr 5545 integer :: nkpt_file 5546 integer :: nblok_d2eig 5547 integer :: iblok_, iblok_d2eig_ 5548 integer :: iband, jband, bandshift, isppol, mband 5549 integer :: ikpt,ipert1,idir1,ipert2,idir2,ii 5550 character(len=500) :: msg 5551 !arrays 5552 integer,allocatable :: flg_d2eig(:,:,:,:) 5553 real(dp) :: qpt(3) 5554 real(dp),allocatable :: nrm(:) 5555 real(dp),allocatable :: matrix_d2eig(:,:,:,:,:,:,:) 5556 real(dp),allocatable :: matrix_d2eig_isppol(:,:,:,:,:,:,:) 5557 5558 ! ************************************************************************ 5559 5560 ncid_d2eig = nctk_idgroup(ncid, 'd2eig') 5561 5562 if (present(iblok_d2eig)) then 5563 iblok_d2eig_= iblok_d2eig 5564 else 5565 ! Recount the blok index 5566 iblok_d2eig_ = 0 5567 do iblok_=1,iblok 5568 if (is_type_d2eig(ddb%typ(iblok_))) then 5569 iblok_d2eig_ = iblok_d2eig_ + 1 5570 end if 5571 end do 5572 end if 5573 5574 ! Sanity check on dimensions 5575 if (MOD(ddb%nband, ddb%nsppol)/=0) then 5576 write(msg,'(a,i5,a,i5)') 'ddb was allocated with nband=',ddb%nband,& 5577 ' but nsppol=',ddb%nsppol 5578 ABI_ERROR(msg) 5579 end if 5580 5581 ! Read kpoints 5582 NCF_CHECK(nctk_get_dim(ncid, "number_of_kpoints", nkpt_file)) 5583 ncerr = nf90_get_var(ncid, nctk_idname(ncid, 'reduced_coordinates_of_kpoints'), ddb%kpt, count=[3,nkpt_file]) 5584 NCF_CHECK(ncerr) 5585 5586 mband = ddb%nband / ddb%nsppol 5587 5588 ncerr = nf90_get_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 5589 'reduced_coordinates_of_qpoints'),& 5590 qpt,& 5591 start=[1,iblok_d2eig_]) 5592 NCF_CHECK(ncerr) 5593 5594 ddb%qpt(:,iblok) = zero 5595 ddb%qpt(1:3,iblok) = qpt 5596 5597 !ncerr = nf90_get_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 5598 ! 'qpoints_normalization'),& 5599 ! nrm(1,iblok),& 5600 ! start=[iblok_d2eig_]) 5601 NCF_CHECK(nctk_get_dim(ncid_d2eig, "number_of_d2eig_blocks", nblok_d2eig)) 5602 ABI_MALLOC(nrm, (nblok_d2eig)) 5603 ncerr = nf90_get_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 5604 'qpoints_normalization'),& 5605 nrm) 5606 NCF_CHECK(ncerr) 5607 ddb%nrm(:,iblok) = zero 5608 ddb%nrm(1,iblok) = nrm(iblok_d2eig_) 5609 ABI_FREE(nrm) 5610 5611 5612 ABI_MALLOC(matrix_d2eig, (2,3,ddb%mpert,3,ddb%mpert,ddb%nband,ddb%nkpt)) 5613 ABI_MALLOC(matrix_d2eig_isppol, (2,3,ddb%mpert,3,ddb%mpert,mband,ddb%nkpt)) 5614 ABI_MALLOC(flg_d2eig, (3,ddb%mpert,3,ddb%mpert)) 5615 5616 do isppol=1,ddb%nsppol 5617 5618 ncerr = nf90_get_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 5619 'matrix_values'),& 5620 matrix_d2eig_isppol,& 5621 start=[1,1,1,1,1,1,1,isppol,iblok_d2eig_]) 5622 NCF_CHECK(ncerr) 5623 5624 if (ddb%nsppol==1) then 5625 matrix_d2eig(:,:,:,:,:,:,:) = matrix_d2eig_isppol(:,:,:,:,:,:,:) 5626 else 5627 bandshift = (isppol - 1) * mband 5628 do ikpt=1,ddb%nkpt 5629 do iband=1,ddb%nband 5630 jband = bandshift + iband 5631 do ipert1=1,ddb%mpert 5632 do idir1=1,3 5633 do ipert2=1,ddb%mpert 5634 do idir2=1,3 5635 do ii=1,2 5636 matrix_d2eig(ii,idir2,ipert2,idir1,ipert1,jband,ikpt)=& 5637 matrix_d2eig_isppol(ii,idir2,ipert2,idir1,ipert1,iband,ikpt) 5638 end do 5639 end do 5640 end do 5641 end do 5642 end do 5643 end do 5644 end do 5645 end if 5646 end do 5647 5648 ncerr = nf90_get_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 5649 'matrix_mask'),& 5650 flg_d2eig,& 5651 start=[1,1,1,1,iblok_d2eig_]) 5652 NCF_CHECK(ncerr) 5653 5654 ! Store values 5655 call ddb%set_d2eig(iblok, matrix_d2eig, flg_d2eig) 5656 5657 ! Free memory 5658 ABI_FREE(matrix_d2eig) 5659 ABI_FREE(matrix_d2eig_isppol) 5660 ABI_FREE(flg_d2eig) 5661 5662 end subroutine ddb_read_d2eig_nc
m_ddb/ddb_read_d2eig_txt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_d2eig_txt
FUNCTION
Read the next DDB block containing second-order derivatives of eigenvalues and store it in block number iblok. The ddb object must have been allocated with nband and nkpt.
INPUTS
unddb=unit for the open file in text format iblok=the block index in the ddb object
OUTPUT
SOURCE
1958 subroutine ddb_read_d2eig_txt(ddb, unddb, iblok) 1959 1960 !Arguments ------------------------------- 1961 !scalars 1962 class(ddb_type),intent(inout) :: ddb 1963 integer, intent(in) :: unddb 1964 integer, intent(in), optional :: iblok 1965 !Local variables ------------------------- 1966 !scalars 1967 integer :: iblok_eig2d 1968 1969 ! ********************************************************************* 1970 1971 iblok_eig2d = 1 1972 if (present(iblok)) iblok_eig2d = iblok 1973 1974 ! GA: Here, nband should really be nband * nsppol. 1975 ! but this is the responsibility of the calling routine 1976 ! see thmeig and merge_ddb 1977 ! FIXME This is inconsistent with ddb_malloc_d2eig... 1978 ! I think I should change this with ddb%nband * ddb%nsppol 1979 call ddb%read_block_txt(iblok_eig2d,ddb%nband,ddb%mpert,ddb%msize,ddb%nkpt,unddb,& 1980 ddb%eig2dval(:,:,:,:),ddb%kpt(:,:)) 1981 1982 end subroutine ddb_read_d2eig_txt
m_ddb/ddb_read_d3E_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_d3E_nc
FUNCTION
Read a DDB block containing 3rd order derivatives of energy.
INPUTS
ncid=netcdf identifier of a file open in reading mode. iblok=index of the block we are setting. iblok_d3E=index of the block we are reading in the d3E group.
OUTPUT
SOURCE
5458 subroutine ddb_read_d3E_nc(ddb, ncid, iblok, iblok_d3E) 5459 5460 !Arguments ------------------------------- 5461 !scalars 5462 class(ddb_type),intent(inout) :: ddb 5463 integer,intent(in) :: ncid,iblok,iblok_d3E 5464 5465 !Local variables ------------------------- 5466 !scalars 5467 integer :: blktyp 5468 integer :: ncid_d3E 5469 integer :: ncerr 5470 !arrays 5471 real(dp) :: qpt(3), nrm(3) 5472 real(dp),allocatable :: matrix_d3E(:,:,:,:,:,:,:) 5473 integer,allocatable :: flg_d3E(:,:,:,:,:,:) 5474 5475 ! ************************************************************************ 5476 5477 ncid_d3E = nctk_idgroup(ncid, 'd3E') 5478 5479 ! Allocate temporary arrays 5480 ABI_MALLOC(matrix_d3E, (2,3,ddb%mpert,3,ddb%mpert,3,ddb%mpert)) 5481 ABI_MALLOC(flg_d3E, (3,ddb%mpert,3,ddb%mpert,3,ddb%mpert)) 5482 5483 ncerr = nf90_get_var(ncid_d3E, nctk_idname(ncid_d3E, 'reduced_coordinates_of_qpoints'),qpt,start=[1,1,iblok_d3E],count=[1,3,1]) 5484 NCF_CHECK(ncerr) 5485 ddb%qpt(1:3,iblok) = qpt(:) 5486 5487 ncerr = nf90_get_var(ncid_d3E, nctk_idname(ncid_d3E, 'reduced_coordinates_of_qpoints'),qpt,start=[2,1,iblok_d3E],count=[1,3,1]) 5488 NCF_CHECK(ncerr) 5489 ddb%qpt(4:6,iblok) = qpt(:) 5490 5491 ncerr = nf90_get_var(ncid_d3E, nctk_idname(ncid_d3E, 'reduced_coordinates_of_qpoints'),qpt,start=[3,1,iblok_d3E],count=[1,3,1]) 5492 NCF_CHECK(ncerr) 5493 ddb%qpt(7:9,iblok) = qpt(:) 5494 5495 ncerr = nf90_get_var(ncid_d3E, nctk_idname(ncid_d3E, 'qpoints_normalization'), nrm, start=[1,iblok_d3E],count=[3,1]) 5496 NCF_CHECK(ncerr) 5497 ddb%nrm(:,iblok) = nrm(:) 5498 5499 NCF_CHECK(nf90_get_var(ncid_d3E, nctk_idname(ncid_d3E, 'matrix_values'), matrix_d3E, start=[1,1,1,1,1,1,1,iblok_d3E])) 5500 NCF_CHECK(nf90_get_var(ncid_d3E, nctk_idname(ncid_d3E, 'matrix_mask'), flg_d3E, start=[1,1,1,1,1,1,iblok_d3E])) 5501 5502 5503 blktyp = ddb%typ(iblok) ! Save block type so it doesnt get overwritten. 5504 5505 call ddb%set_d3matr(iblok, matrix_d3E, flg_d3E) 5506 5507 ddb%typ(iblok) = blktyp 5508 5509 ! Free memory 5510 ABI_FREE(matrix_d3E) 5511 ABI_FREE(flg_d3E) 5512 5513 end subroutine ddb_read_d3E_nc
m_ddb/ddb_read_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_nc
FUNCTION
This subroutine reads data from the DDB.nc file and constructs an instance of ddb_type It also returns an instance of crystal_t with the crystalline structure reported in the DDB file and the DDB header.
INPUTS
filename=DDB filename. comm=MPI communicator. prtvol=Verbosity level raw = 1 -> do not perform any symetrization or transformation to cartesian coordinates. 0 (default) -> do perform these transformations.
OUTPUT
ddb<type(ddb_type)>=Object storing the DDB results. crystal<type(crystal_t)>=Crystal structure parameters ddb_hdr= Header of the DDB file.
SOURCE
2659 subroutine ddb_read_nc(ddb, filename, ddb_hdr, crystal, comm, prtvol, raw) 2660 2661 !Arguments ------------------------------- 2662 !scalars 2663 class(ddb_type),intent(inout) :: ddb 2664 integer,intent(in) :: comm 2665 integer,optional,intent(in) :: prtvol, raw 2666 character(len=*),intent(in) :: filename 2667 type(crystal_t),intent(out) :: crystal 2668 type(ddb_hdr_type),intent(out) :: ddb_hdr 2669 !array 2670 2671 !Local variables------------------------------- 2672 !scalars 2673 integer,parameter :: master=0 2674 integer :: prtvol_, raw_ 2675 integer :: ncid 2676 integer :: iblok,iblok_d0E,iblok_d1E,iblok_d2E,iblok_d3E,iblok_d2eig 2677 2678 !arrays 2679 !character(len=132),allocatable :: title(:) 2680 2681 ! ************************************************************************ 2682 2683 DBG_ENTER("COLL") 2684 2685 if (present(raw)) then 2686 raw_ = raw 2687 else 2688 raw_ = 0 2689 end if 2690 2691 ! GA: Not really used so far 2692 if (present(prtvol)) then 2693 prtvol_ = prtvol 2694 else 2695 prtvol_ = 0 2696 end if 2697 2698 ! Read header 2699 call ddb_hdr%open_read_nc(filename, comm) 2700 ncid = ddb_hdr%ncid 2701 2702 if (xmpi_comm_rank(comm) == master) then 2703 2704 ddb%nsppol = ddb_hdr%nsppol 2705 2706 ! Copy dimensions from header and allocate arrays 2707 call ddb%malloc(ddb_hdr%msize, ddb_hdr%nblok, ddb_hdr%natom, & 2708 ddb_hdr%ntypat, ddb_hdr%mpert,& 2709 ddb_hdr%nkpt, ddb_hdr%mband) 2710 2711 ! Copy arrays from header 2712 ddb%typ(:) = ddb_hdr%typ(:) 2713 ddb%amu(:) = ddb_hdr%crystal%amu(:) 2714 ddb%acell(:) = one 2715 ddb%rprim(:,:) = ddb_hdr%crystal%rprimd(:,:) 2716 ddb%gprim(:,:) = ddb_hdr%crystal%gprimd(:,:) 2717 2718 ! --------------- 2719 ! Read all blocks 2720 ! --------------- 2721 iblok_d0E = 0 2722 iblok_d1E = 0 2723 iblok_d2E = 0 2724 iblok_d3E = 0 2725 iblok_d2eig = 0 2726 2727 do iblok=1,ddb%nblok 2728 2729 if (is_type_d0E(ddb%typ(iblok))) then 2730 iblok_d0E = iblok_d0E + 1 2731 call ddb%read_d0E_nc(ncid, iblok, iblok_d0E) 2732 2733 else if (is_type_d1E(ddb%typ(iblok))) then 2734 iblok_d1E = iblok_d1E + 1 2735 call ddb%read_d1E_nc(ncid, iblok, iblok_d1E) 2736 2737 else if (is_type_d2E(ddb%typ(iblok))) then 2738 iblok_d2E = iblok_d2E + 1 2739 call ddb%read_d2E_nc(ncid, iblok, iblok_d2E) 2740 2741 else if (is_type_d3E(ddb%typ(iblok))) then 2742 iblok_d3E = iblok_d3E + 1 2743 call ddb%read_d3E_nc(ncid, iblok, iblok_d3E) 2744 2745 else if (is_type_d2eig(ddb%typ(iblok))) then 2746 iblok_d2eig = iblok_d2eig + 1 2747 ! GA: It is kind of weird to call this function inside a loop, 2748 ! because the ddb can only hold a single block of d2eig data. 2749 call ddb%read_d2eig_nc(ncid, iblok, iblok_d2eig) 2750 2751 end if 2752 2753 ! Symmetrize and transform if raw==0 2754 if (raw_==0) then 2755 call ddb%symmetrize_and_transform(ddb_hdr%crystal,iblok) 2756 end if 2757 2758 end do 2759 2760 end if 2761 2762 ! Close the file 2763 call ddb_hdr%close() 2764 2765 ! -------------- 2766 ! Broadcast data 2767 ! -------------- 2768 if (xmpi_comm_size(comm) > 1) then 2769 call ddb%bcast(comm) 2770 call ddb_hdr%bcast(comm) 2771 end if 2772 2773 ! Copy crystal 2774 call ddb_hdr%crystal%copy(crystal) 2775 2776 DBG_EXIT("COLL") 2777 2778 end subroutine ddb_read_nc
m_ddb/ddb_read_txt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_read_txt
FUNCTION
This subroutine reads data from the DDB file and constructs an instance of ddb_type It also returns an instance of crystal_t with the crystalline structure reported in the DDB file and the DDB header.
INPUTS
filename=DDB filename. comm=MPI communicator. prtvol=Verbosity level raw = 1 -> do not perform any symetrization or transformation to cartesian coordinates. 0 (default) -> do perform these transformations.
OUTPUT
ddb<type(ddb_type)>=Object storing the DDB results. crystal<type(crystal_t)>=Crystal structure parameters ddb_hdr= Header of the DDB file.
SOURCE
2473 subroutine ddb_read_txt(ddb, filename, ddb_hdr, crystal, comm, prtvol, raw) 2474 2475 !Arguments ------------------------------- 2476 !scalars 2477 class(ddb_type),intent(inout) :: ddb 2478 integer,intent(in) :: comm 2479 integer,optional,intent(in) :: prtvol, raw 2480 character(len=*),intent(in) :: filename 2481 type(crystal_t),intent(out) :: Crystal 2482 type(ddb_hdr_type),intent(out) :: ddb_hdr 2483 2484 !Local variables------------------------------- 2485 !scalars 2486 integer,parameter :: master=0 2487 integer :: msym,dimekb,lmnmax,mband,nkpt,ntypat,nsym,usepaw 2488 integer :: mpert,msize,natom,nblok,occopt,nsppol 2489 real(dp) :: ucvol 2490 2491 !arrays 2492 integer,allocatable :: symrec(:,:,:),symrel(:,:,:),symafm(:),indsym(:,:,:),typat(:) 2493 real(dp) :: acell(3),gmet(3,3),gprim(3,3),rmet(3,3),rprim(3,3) 2494 real(dp),allocatable :: amu(:),xcart(:),xred(:,:),zion(:),znucl(:),tnons(:,:) 2495 2496 ! ************************************************************************ 2497 2498 DBG_ENTER("COLL") 2499 2500 ! Must read natom from the DDB before being able to allocate some arrays needed for invars9 2501 call ddb_hdr%open_read_txt(filename, comm) 2502 2503 ! GA: clean this up. Not all of it is useful 2504 nblok = ddb_hdr%nblok 2505 msym = ddb_hdr%msym 2506 natom = ddb_hdr%natom 2507 ntypat = ddb_hdr%ntypat 2508 mband = ddb_hdr%mband 2509 nkpt = ddb_hdr%nkpt 2510 nsppol = ddb_hdr%nsppol 2511 usepaw = ddb_hdr%usepaw 2512 dimekb = ddb_hdr%psps%dimekb 2513 lmnmax = ddb_hdr%psps%lmnmax 2514 occopt = ddb_hdr%occopt 2515 mpert = ddb_hdr%mpert 2516 msize = ddb_hdr%msize 2517 2518 ! Master reads and then broadcasts data. 2519 if (xmpi_comm_rank(comm) == master) then 2520 2521 ! Allocate arrays depending on msym (which is actually fixed to nsym inside inprep8) 2522 ABI_MALLOC(symrel,(3,3,msym)) 2523 ABI_MALLOC(symafm,(msym)) 2524 ABI_MALLOC(tnons,(3,msym)) 2525 ABI_MALLOC(typat,(natom)) 2526 ABI_MALLOC(xred,(3,natom)) 2527 ABI_MALLOC(zion,(ntypat)) 2528 ABI_MALLOC(znucl,(ntypat)) 2529 2530 ABI_MALLOC(symrec,(3,3,msym)) 2531 ABI_MALLOC(indsym,(4,msym,natom)) 2532 ABI_MALLOC(xcart,(3*natom)) 2533 ABI_MALLOC(amu,(ntypat)) 2534 2535 ddb%nsppol = nsppol 2536 call ddb%malloc(msize, nblok, natom, ntypat, mpert) 2537 2538 ! GA: FIXME 2539 ! Should clean this up. Lots of the arguments are not needed. 2540 ! In particular, rprim and acell could be taken from ddb_hdr%crystal 2541 ! which is already initialized at this point 2542 call rdddb9(ddb, ddb_hdr, ddb_hdr%unddb,& 2543 acell,amu,gmet,gprim,indsym,& 2544 mband,mpert,msize,msym,& 2545 natom,nkpt,nsym,ntypat,& 2546 rmet,rprim,symrec,symrel,symafm,& 2547 tnons,typat,ucvol,xcart,xred,zion,znucl,raw) 2548 2549 ABI_FREE(symrec) 2550 ABI_FREE(indsym) 2551 ABI_FREE(xcart) 2552 2553 ! Save variables needed to call legacy code. 2554 ddb%acell = acell 2555 ddb%rprim = rprim 2556 ddb%gprim = gprim 2557 2558 !call ddb%set_brav(brav) 2559 2560 ! Other useful quantities. 2561 ! 2 is to preserve the old behaviour 2562 ddb%prtvol = 2; if (present(prtvol)) ddb%prtvol = prtvol 2563 ddb%occopt = occopt 2564 ddb%amu = amu 2565 ABI_FREE(amu) 2566 2567 ! These were not needed because crystal is already initialized in the header. 2568 ABI_FREE(symrel) 2569 ABI_FREE(symafm) 2570 ABI_FREE(tnons) 2571 ABI_FREE(typat) 2572 ABI_FREE(xred) 2573 ABI_FREE(zion) 2574 ABI_FREE(znucl) 2575 2576 end if 2577 2578 call ddb_hdr%close() 2579 2580 if (xmpi_comm_size(comm) > 1) then 2581 call ddb%bcast(comm) 2582 call ddb_hdr%bcast(comm) 2583 2584 !! GA: This seems superfluous now... 2585 !call xmpi_bcast(nsym, master, comm, ierr) 2586 !call xmpi_bcast(symrel, master, comm, ierr) 2587 !call xmpi_bcast(symafm, master, comm, ierr) 2588 !call xmpi_bcast(typat, master, comm, ierr) 2589 !call xmpi_bcast(acell, master, comm, ierr) 2590 !call xmpi_bcast(occopt, master, comm, ierr) 2591 !call xmpi_bcast(gprim, master, comm, ierr) 2592 !call xmpi_bcast(rprim, master, comm, ierr) 2593 !call xmpi_bcast(tnons, master, comm, ierr) 2594 !call xmpi_bcast(xred, master, comm, ierr) 2595 !call xmpi_bcast(zion, master, comm, ierr) 2596 !call xmpi_bcast(znucl, master, comm, ierr) 2597 end if 2598 2599 call ddb_hdr%crystal%copy(Crystal) 2600 2601 !! Initialize crystal_t object. 2602 !call mkrdim(acell,rprim,rprimd) 2603 2604 !! GA: These variables are hardcoded which means the crystal object 2605 !! is not reliable for antiferro systems or alchemical potentials 2606 !! when it is read from a text DDB file. 2607 !npsp = ntypat; space_group = 0; timrev = 2 2608 !use_antiferro=.FALSE. !; use_antiferro=(nspden==2.and.nsppol==1) 2609 !ABI_MALLOC(title, (ntypat)) 2610 2611 !do ii=1,ntypat 2612 ! write(title(ii),'(a,i0)')"No title for typat ",ii 2613 !end do 2614 2615 !! Warning znucl is dimensioned with ntypat = nspsp hence alchemy is not supported here 2616 !call crystal_init(ddb%amu,Crystal,space_group,natom,npsp,ntypat,nsym,rprimd,typat,xred,& 2617 ! zion,znucl,timrev,use_antiferro,.FALSE.,title,& 2618 ! symrel=symrel(:,:,1:nsym),tnons=tnons(:,1:nsym),symafm=symafm(1:nsym)) 2619 2620 !ABI_FREE(title) 2621 !ABI_FREE(symrel) 2622 !ABI_FREE(symafm) 2623 !ABI_FREE(tnons) 2624 !ABI_FREE(typat) 2625 !ABI_FREE(xred) 2626 !ABI_FREE(zion) 2627 !ABI_FREE(znucl) 2628 2629 DBG_EXIT("COLL") 2630 2631 end subroutine ddb_read_txt
m_ddb/ddb_set_brav [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_brav
FUNCTION
Modify the current values of rprim according to bravais lattice. Perform some checks on the primitive vectors before rescaling them such that rprim(1,2)=0.5
INPUTS
brav
1 -> No rescaling.
other -> Check and rescale.
Note that the meaning of brav is
1 or -1 -> simple lattice
2 -> face-centered cubic
3 -> body-centered lattice
4 -> hexagonal lattice (D6h)
OUTPUT
SOURCE
1166 subroutine ddb_set_brav(ddb, brav) 1167 1168 !Arguments ------------------------------- 1169 !array 1170 class(ddb_type),intent(inout) :: ddb 1171 !scalars 1172 integer,intent(in) :: brav 1173 1174 !Local variables------------------------------- 1175 !scalars 1176 real(dp) :: factor 1177 character(len=500) :: msg 1178 1179 ! ************************************************************************* 1180 1181 ! Renormalize rprim to possibly satisfy the constraint abs(rprim(1,2))=half when abs(brav)/=1 1182 ! This section is needed to preserver the behaviour of the old implementation. 1183 if (abs(brav)/=1 .and. abs(abs(ddb%rprim(1,2))-half)>tol10) then 1184 if(abs(ddb%rprim(1,2))<tol6)then 1185 write(msg, '(a,i0,7a)' )& 1186 'The input DDB value of brav is ',brav,',',ch10,& 1187 'and the one of rprim(1,2) is zero.',ch10,& 1188 'These are incompatible',ch10,& 1189 'Action: check the value of brav and rprim(1,2) in your DDB.' 1190 ABI_ERROR(msg) 1191 end if 1192 factor = abs(ddb%rprim(1,2)) * two 1193 ddb%acell(:) = ddb%acell(:) * factor 1194 ddb%rprim(:,:) = ddb%rprim(:,:) / factor 1195 ddb%gprim(:,:) = ddb%gprim(:,:) * factor 1196 end if 1197 1198 end subroutine ddb_set_brav
m_ddb/ddb_set_d1matr [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_d1matr
FUNCTION
Set values for the first-order derivative matrix.
INPUTS
iblok=index of the block being set. d1matr=the first-order derivative matrix. flg=flag to indicate presence of a given element.
SOURCE
1075 subroutine ddb_set_d1matr(ddb, iblok, d1matr, flg) 1076 1077 !Arguments ------------------------------- 1078 !array 1079 class(ddb_type),intent(inout) :: ddb 1080 real(dp), intent(in) :: d1matr(2,3,ddb%mpert) 1081 integer, intent(in) :: flg(3,ddb%mpert) 1082 !scalars 1083 integer,intent(in) :: iblok 1084 1085 !Local variables ------------------------- 1086 !scalars 1087 integer :: ii,ipert1,idir1 1088 1089 ! ************************************************************************ 1090 1091 ii=0 1092 do ipert1=1,ddb%mpert 1093 do idir1=1,3 1094 ii=ii+1 1095 ddb%val(1,ii,iblok) = d1matr(1,idir1,ipert1) 1096 ddb%val(2,ii,iblok) = d1matr(2,idir1,ipert1) 1097 ddb%flg(ii,iblok) = flg(idir1,ipert1) 1098 end do 1099 end do 1100 1101 end subroutine ddb_set_d1matr
m_ddb/ddb_set_d2eig [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_d2eig
FUNCTION
Set values for the second-order derivatives of eigenvalues.
INPUTS
iblok=index of the block we are setting. d2eig=the second-order derivative of eigenvalues. flg=flag to indicate presence of a given element.
OUTPUT
NOTE Does not handle spin index. Also, sometimes, d2eig is available with flat index
SOURCE
5884 subroutine ddb_set_d2eig(ddb, iblok, d2eig, flg) 5885 5886 !Arguments ------------------------------- 5887 !array 5888 class(ddb_type),intent(inout) :: ddb 5889 real(dp),intent(in) :: d2eig(2,3,ddb%mpert,3,ddb%mpert,ddb%nband,ddb%nkpt) 5890 integer,intent(in) :: flg(3,ddb%mpert,3,ddb%mpert) 5891 !scalars 5892 integer,intent(in) :: iblok 5893 5894 !Local variables ------------------------- 5895 !scalars 5896 integer :: idir1,idir2,ipert1,ipert2,index,iband,ikpt,ii 5897 5898 ! ************************************************************************ 5899 5900 do ikpt=1,ddb%nkpt 5901 do iband=1,ddb%nband 5902 ii = 0 5903 do ipert2=1,ddb%mpert 5904 do idir2=1,3 5905 do ipert1=1,ddb%mpert 5906 do idir1=1,3 5907 index=idir1+3*((ipert1-1)+ddb%mpert*((idir2-1)+3*(ipert2-1))) 5908 ddb%eig2dval(1,index,iband,ikpt)=d2eig(1,idir1,ipert1,idir2,ipert2,iband,ikpt) 5909 ddb%eig2dval(2,index,iband,ikpt)=d2eig(2,idir1,ipert1,idir2,ipert2,iband,ikpt) 5910 ddb%flg(index,iblok)=flg(idir1,ipert1,idir2,ipert2) 5911 end do !idir1 5912 end do !pert1 5913 end do !idir2 5914 end do !pert2 5915 end do !band 5916 end do !kpt 5917 5918 end subroutine ddb_set_d2eig
m_ddb/ddb_set_d2eig_reshape [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_d2eig_reshape
FUNCTION
Set values for the second-order derivatives of eigenvalues and reshape the array received.
INPUTS
iblok=index of the block we are setting. d2eig=the second-order derivative of eigenvalues. flg=flag to indicate presence of a given element. blktyp=block type 5->real part 6->imaginary part (broadening)
OUTPUT
SOURCE
5943 subroutine ddb_set_d2eig_reshape(ddb, iblok, d2eig, flg, blktyp) 5944 5945 !Arguments ------------------------------- 5946 !array 5947 class(ddb_type),intent(inout) :: ddb 5948 real(dp),intent(in) :: d2eig(2,ddb%nband*ddb%nsppol,ddb%nkpt,3,ddb%mpert,3,ddb%mpert) 5949 integer,intent(in) :: flg(3,ddb%mpert,3,ddb%mpert) 5950 !scalars 5951 integer,intent(in) :: iblok 5952 integer,intent(in),optional :: blktyp 5953 5954 !Local variables ------------------------- 5955 !scalars 5956 integer :: idir1,idir2,ipert1,ipert2,index,iband,ikpt,ii,mband 5957 5958 ! ************************************************************************ 5959 5960 ddb%typ(iblok) = BLKTYP_d2eig_re 5961 if (present(blktyp)) ddb%typ(iblok) = blktyp 5962 5963 ! Spin polarization is wrapped in band number 5964 mband = ddb%nband * ddb%nsppol 5965 5966 do ikpt=1,ddb%nkpt 5967 do iband=1,mband 5968 ii = 0 5969 do ipert2=1,ddb%mpert 5970 do idir2=1,3 5971 do ipert1=1,ddb%mpert 5972 do idir1=1,3 5973 index=idir1+3*((ipert1-1)+ddb%mpert*((idir2-1)+3*(ipert2-1))) 5974 ddb%flg(index,iblok)=flg(idir1,ipert1,idir2,ipert2) 5975 if (ddb%flg(index,iblok) > 0) then 5976 ddb%eig2dval(1,index,iband,ikpt)=d2eig(1,iband,ikpt,idir1,ipert1,idir2,ipert2) 5977 ddb%eig2dval(2,index,iband,ikpt)=d2eig(2,iband,ikpt,idir1,ipert1,idir2,ipert2) 5978 end if 5979 end do !idir1 5980 end do !pert1 5981 end do !idir2 5982 end do !pert2 5983 end do !band 5984 end do !kpt 5985 5986 end subroutine ddb_set_d2eig_reshape
m_ddb/ddb_set_d2matr [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_d2matr
FUNCTION
Set values for the second-order derivative matrix.
INPUTS
iblok=index of the block being set. d2matr=the second-order derivative matrix. flg=flag to indicate presence of a given element.
OUTPUT
SOURCE
777 subroutine ddb_set_d2matr(ddb, iblok, d2matr, flg) 778 779 !Arguments ------------------------------- 780 class(ddb_type),intent(inout) :: ddb 781 !scalars 782 integer,intent(in) :: iblok 783 !array 784 real(dp), intent(in) :: d2matr(2,3,ddb%mpert,3,ddb%mpert) 785 integer, intent(in) :: flg(3,ddb%mpert,3,ddb%mpert) 786 787 !Local variables ------------------------- 788 !scalars 789 integer :: idir1,idir2,ii,ipert1,ipert2 790 791 ! ************************************************************************ 792 793 ii=0 794 do ipert2=1,ddb%mpert 795 do idir2=1,3 796 do ipert1=1,ddb%mpert 797 do idir1=1,3 798 ii=ii+1 799 ddb%flg(ii,iblok) = flg(idir1,ipert1,idir2,ipert2) 800 ddb%val(1,ii,iblok) = d2matr(1,idir1,ipert1,idir2,ipert2) 801 ddb%val(2,ii,iblok) = d2matr(2,idir1,ipert1,idir2,ipert2) 802 end do 803 end do 804 end do 805 end do 806 807 end subroutine ddb_set_d2matr
m_ddb/ddb_set_d3matr [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_d3matr
FUNCTION
Set values for the third-order derivative matrix.
INPUTS
iblok=index of the block we are setting. d2matr=the third-order derivative matrix. flg=flag to indicate presence of a given element. lw=whether this type of perturbation correspond to longwave derivatives
OUTPUT
SOURCE
5684 subroutine ddb_set_d3matr(ddb, iblok, d3matr, flg, lw) 5685 5686 !Arguments ------------------------------- 5687 !array 5688 class(ddb_type),intent(inout) :: ddb 5689 real(dp),intent(in) :: d3matr(2,3,ddb%mpert,3,ddb%mpert,3,ddb%mpert) 5690 integer,intent(in) :: flg(3,ddb%mpert,3,ddb%mpert,3,ddb%mpert) 5691 !scalars 5692 integer,intent(in) :: iblok 5693 logical,intent(in),optional :: lw 5694 5695 !Local variables ------------------------- 5696 !scalars 5697 integer :: idir1,idir2,idir3,ipert1,ipert2,ipert3,index,mpert 5698 5699 ! ************************************************************************ 5700 5701 mpert = ddb%mpert 5702 5703 ! GA: Should be consistent among all ddb_set_dxmatr routines 5704 ! and alway have options to specify block type... 5705 ddb%typ(iblok) = BLKTYP_d3E_xx 5706 if (present(lw)) then 5707 if (lw) then 5708 ddb%typ(iblok) = BLKTYP_d3E_lw 5709 end if 5710 end if 5711 5712 do ipert1=1,mpert 5713 do idir1=1,3 5714 do ipert2=1,mpert 5715 do idir2=1,3 5716 do ipert3=1,mpert 5717 do idir3=1,3 5718 5719 index = idir1 + 3*((ipert1-1)+mpert*((idir2-1) + & 5720 & 3*((ipert2-1)+mpert*((idir3-1) + 3*(ipert3-1))))) 5721 5722 ddb%flg(index,iblok) = flg(idir1,ipert1,idir2,ipert2,idir3,ipert3) 5723 ddb%val(:,index,iblok)= d3matr(:,idir1,ipert1,idir2,ipert2,idir3,ipert3) 5724 5725 end do 5726 end do 5727 end do 5728 end do 5729 end do 5730 end do 5731 5732 5733 end subroutine ddb_set_d3matr
m_ddb/ddb_set_etotal [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_etotal
FUNCTION
Set the total energy
INPUTS
etotal=the total energy. iblok=index of the block we are setting.
OUTPUT
SOURCE
1121 subroutine ddb_set_etotal(ddb, etotal, iblok) 1122 1123 !Arguments ------------------------------- 1124 !array 1125 class(ddb_type),intent(inout) :: ddb 1126 !scalars 1127 real(dp), intent(in) :: etotal 1128 integer,intent(in) :: iblok 1129 1130 ! ************************************************************************ 1131 1132 ddb%typ(iblok) = BLKTYP_d0E_xx 1133 ddb%val(1,1,iblok) = etotal 1134 ddb%val(2,1,iblok) = zero 1135 ddb%flg(1,iblok) = 1 1136 1137 end subroutine ddb_set_etotal
m_ddb/ddb_set_gred [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_gred
FUNCTION
Set the forces in reduced coordinates (Hartree).
INPUTS
gred=the gradient of the total energy with respect to change of reduced coordinates iblok=index of the block being set.
OUTPUT
SOURCE
887 subroutine ddb_set_gred(ddb, gred, iblok) 888 889 !Arguments ------------------------------- 890 !array 891 class(ddb_type),intent(inout) :: ddb 892 real(dp), intent(in) :: gred(3,ddb%natom) 893 !scalars 894 integer,intent(in) :: iblok 895 896 !Local variables ------------------------- 897 !scalars 898 integer :: idir, iatom, indx 899 900 ! ************************************************************************ 901 902 ddb%typ(iblok) = BLKTYP_d1E_xx 903 indx = 0 904 do iatom = 1, ddb%natom 905 do idir = 1, 3 906 indx = indx + 1 907 ddb%flg(indx,iblok) = 1 908 ddb%val(1,indx,iblok) = gred(idir,iatom) 909 ddb%val(2,indx,iblok) = zero 910 end do 911 end do 912 913 end subroutine ddb_set_gred
m_ddb/ddb_set_pel [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_gred
FUNCTION
Set the electronic polarization.
INPUTS
pel=ucvol times the electronic polarization in reduced coordinates. flg=flag to indicate presence of a given element. iblok=index of the block being set.
OUTPUT
SOURCE
934 subroutine ddb_set_pel(ddb, pel, flg, iblok) 935 936 !Arguments ------------------------------- 937 !array 938 class(ddb_type),intent(inout) :: ddb 939 real(dp), intent(in) :: pel(3) 940 integer,intent(in) :: flg(3) 941 !scalars 942 integer,intent(in) :: iblok 943 944 !Local variables ------------------------- 945 !scalars 946 integer :: idir, indx 947 948 ! ************************************************************************ 949 950 ddb%typ(iblok) = BLKTYP_d1E_xx 951 indx = 3*ddb%natom + 3 952 do idir = 1, 3 953 indx = indx + 1 954 ddb%flg(indx,iblok) = flg(idir) 955 ddb%val(1,indx,iblok) = pel(idir) 956 ddb%val(2,indx,iblok) = zero 957 end do 958 959 end subroutine ddb_set_pel
m_ddb/ddb_set_qpt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_qpt
FUNCTION
Set the q-point wavevector for a certain block. In case of 3rd order derivatives, three q-points need to be specified with the constrain q1 + q2 + q3 = 0.
INPUTS
iblok=index of the block being set. qpt=reduced coordinates of first qpoint qpt2=reduced coordinates of second qpoint qpt3=reduced coordinates of third qpoint
OUTPUT
SOURCE
730 subroutine ddb_set_qpt(ddb, iblok, qpt, qpt2, qpt3) 731 732 !Arguments ------------------------------- 733 !array 734 class(ddb_type),intent(inout) :: ddb 735 real(dp), intent(in) :: qpt(3) 736 real(dp), intent(in),optional :: qpt2(3) 737 real(dp), intent(in),optional :: qpt3(3) 738 !scalars 739 integer,intent(in) :: iblok 740 741 ! ************************************************************************ 742 743 ddb%qpt(1:3,iblok) = qpt(1:3) 744 ddb%nrm(1,iblok) = one 745 746 if (present(qpt2)) then 747 ddb%qpt(4:6,iblok) = qpt2(1:3) 748 ddb%nrm(2,iblok) = one 749 end if 750 751 if (present(qpt3)) then 752 ddb%qpt(7:9,iblok) = qpt3(1:3) 753 ddb%nrm(3,iblok) = one 754 end if 755 756 end subroutine ddb_set_qpt
m_ddb/ddb_set_strten [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_strten
FUNCTION
Set the stress tensor.
INPUTS
strten=the stress tensor in cartesian coordinates. iblok=index of the block we are setting.
OUTPUT
SOURCE
979 subroutine ddb_set_strten(ddb, strten, iblok) 980 981 !Arguments ------------------------------- 982 !array 983 class(ddb_type),intent(inout) :: ddb 984 real(dp), intent(in) :: strten(6) 985 !scalars 986 integer,intent(in) :: iblok 987 988 !Local variables ------------------------- 989 !scalars 990 integer :: indx 991 992 ! ************************************************************************ 993 994 ddb%typ(iblok) = BLKTYP_d1E_xx 995 indx = 3*ddb%natom + 6 996 997 ddb%flg(indx+1:indx+6,1) = 1 998 ddb%val(1,indx+1:indx+6,1) = strten(1:6) 999 ddb%val(2,indx+1:indx+6,1) = zero 1000 1001 end subroutine ddb_set_strten
m_ddb/ddb_set_typ [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_set_typ
FUNCTION
Set the blok typ for one block
INPUTS
iblok: block index typ: type of block
OUTPUT
SOURCE
1218 subroutine ddb_set_typ(ddb, iblok, typ) 1219 1220 !Arguments ------------------------------- 1221 !array 1222 class(ddb_type),intent(inout) :: ddb 1223 !scalars 1224 integer,intent(in) :: iblok,typ 1225 ! ************************************************************************* 1226 1227 ddb%typ(iblok) = typ 1228 1229 end subroutine ddb_set_typ
m_ddb/ddb_symmetrize_and_transform [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_symmetrize_and_transform
FUNCTION
First apply symmetry operations, then transform the second-derivative matrix from reduced coordinates to cartesian coordinates, and also 1) add the ionic part of the effective charges, 2) normalize the electronic dielectric tensor, and add the vacuum polarisation
INPUTS
crystal<type(crystal_t)> iblock=the block index on which to act
SIDE EFFECTS
ddb<type(ddb_type)>=
OUTPUT
SOURCE
4212 subroutine ddb_symmetrize_and_transform(ddb, crystal, iblok) 4213 4214 !Arguments ------------------------------- 4215 !scalars 4216 class(ddb_type),intent(inout) :: ddb 4217 class(crystal_t),intent(in) :: crystal 4218 integer,intent(in) :: iblok 4219 !integer,intent(inout) :: indsym(4,msym,natom) 4220 !integer,intent(out) :: symrec(3,3,msym),symrel(3,3,msym),symafm(msym) 4221 4222 !Local variables------------------------------- 4223 !scalars 4224 integer :: mpert,natom,nsym,ntypat 4225 integer :: nsize,timrev 4226 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 4227 !arrays 4228 !integer :: symq(4,2,msym) 4229 integer,allocatable :: symq(:,:,:) 4230 integer,allocatable :: car3flg(:,:,:,:,:,:),carflg(:,:,:,:) 4231 integer,allocatable :: tmpflg(:,:,:,:,:,:),rfpert(:,:,:,:,:,:) 4232 real(dp) :: gprimd(3,3),qpt(3),rprimd(3,3) 4233 real(dp),allocatable :: d2cart(:,:,:,:,:),d3cart(:,:,:,:,:,:,:) 4234 real(dp),allocatable :: tmpval(:,:,:,:,:,:,:) 4235 4236 ! ************************************************************************ 4237 4238 mpert = ddb%mpert 4239 natom = ddb%natom 4240 ntypat = crystal%ntypat 4241 4242 nsym = crystal%nsym 4243 rprimd(:,:) = crystal%rprimd(:,:) 4244 gprimd(:,:) = crystal%gprimd(:,:) 4245 4246 ABI_MALLOC(symq,(4,2,nsym)) 4247 4248 ! Here complete the matrix by symmetrisation of the existing elements 4249 if (ddb%typ(iblok) == BLKTYP_d2E_ns .or. ddb%typ(iblok) == BLKTYP_d2E_st) then 4250 4251 qpt(1)=ddb%qpt(1,iblok)/ddb%nrm(1,iblok) 4252 qpt(2)=ddb%qpt(2,iblok)/ddb%nrm(1,iblok) 4253 qpt(3)=ddb%qpt(3,iblok)/ddb%nrm(1,iblok) 4254 4255 ! Examine the symmetries of the q wavevector 4256 call littlegroup_q(crystal%nsym,qpt,symq,crystal%symrec,crystal%symafm,timrev,prtvol=0) 4257 4258 !GA: Note that d2sym3 and cart29 expect different shapes for tmpflg and tmpval 4259 ! hence the extra dimensions 4260 nsize=3*mpert*3*mpert 4261 ABI_MALLOC(tmpflg,(3,mpert,3,mpert,1,1)) 4262 ABI_MALLOC(tmpval,(2,3,mpert,3,mpert,1,1)) 4263 4264 tmpflg(:,:,:,:,1,1) = reshape(ddb%flg(1:nsize,iblok), shape = (/3,mpert,3,mpert/)) 4265 tmpval(1,:,:,:,:,1,1) = reshape(ddb%val(1,1:nsize,iblok), shape = (/3,mpert,3,mpert/)) 4266 tmpval(2,:,:,:,:,1,1) = reshape(ddb%val(2,1:nsize,iblok), shape = (/3,mpert,3,mpert/)) 4267 4268 ! Then apply symmetry operations 4269 call d2sym3(tmpflg,tmpval,crystal%indsym,mpert,natom,nsym,qpt,symq,crystal%symrec,crystal%symrel,timrev,1) 4270 4271 ! Transform the dynamical matrix in cartesian coordinates 4272 ABI_MALLOC(carflg,(3,mpert,3,mpert)) 4273 ABI_MALLOC(d2cart,(2,3,mpert,3,mpert)) 4274 4275 call cart29(tmpflg,tmpval,carflg,d2cart,gprimd,1,mpert,natom,1,ntypat,rprimd,crystal%typat,crystal%ucvol,crystal%zion) 4276 4277 ddb%flg(1:nsize,iblok) = reshape(carflg,shape = (/3*mpert*3*mpert/)) 4278 ddb%val(1,1:nsize,iblok) = reshape(d2cart(1,:,:,:,:), shape = (/3*mpert*3*mpert/)) 4279 ddb%val(2,1:nsize,iblok) = reshape(d2cart(2,:,:,:,:), shape = (/3*mpert*3*mpert/)) 4280 4281 ABI_FREE(carflg) 4282 ABI_FREE(d2cart) 4283 ABI_FREE(tmpflg) 4284 ABI_FREE(tmpval) 4285 4286 else if (ddb%typ(iblok) == BLKTYP_d3E_xx) then 4287 4288 nsize=3*mpert*3*mpert*3*mpert 4289 ABI_MALLOC(tmpflg,(3,mpert,3,mpert,3,mpert)) 4290 ABI_MALLOC(tmpval,(2,3,mpert,3,mpert,3,mpert)) 4291 ABI_MALLOC(rfpert,(3,mpert,3,mpert,3,mpert)) 4292 4293 tmpflg(:,:,:,:,:,:) = reshape(ddb%flg(1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/)) 4294 tmpval(1,:,:,:,:,:,:) = reshape(ddb%val(1,1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/)) 4295 tmpval(2,:,:,:,:,:,:) = reshape(ddb%val(2,1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/)) 4296 4297 ! Set the elements that are zero by symmetry for raman and 4298 ! non-linear optical susceptibility tensors 4299 rfpert = 0 4300 rfpert(:,natom+2,:,natom+2,:,natom+2) = 1 4301 rfpert(:,1:natom,:,natom+2,:,natom+2) = 1 4302 rfpert(:,natom+2,:,1:natom,:,natom+2) = 1 4303 rfpert(:,natom+2,:,natom+2,:,1:natom) = 1 4304 call sytens(crystal%indsym,mpert,natom,nsym,rfpert,crystal%symrec,crystal%symrel) 4305 do i1pert = 1,mpert 4306 do i2pert = 1,mpert 4307 do i3pert = 1,mpert 4308 do i1dir=1,3 4309 do i2dir=1,3 4310 do i3dir=1,3 4311 if ((rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) .and. & 4312 (tmpflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)) then 4313 tmpval(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = zero 4314 tmpflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)=1 4315 end if 4316 end do 4317 end do 4318 end do 4319 end do 4320 end do 4321 end do 4322 4323 call d3sym(tmpflg,tmpval,crystal%indsym,mpert,natom,nsym,crystal%symrec,crystal%symrel) 4324 4325 ABI_MALLOC(d3cart,(2,3,mpert,3,mpert,3,mpert)) 4326 ABI_MALLOC(car3flg,(3,mpert,3,mpert,3,mpert)) 4327 4328 call nlopt(tmpflg,car3flg,tmpval,d3cart,gprimd,mpert,natom,rprimd,crystal%ucvol) 4329 4330 ddb%flg(1:nsize,iblok) = reshape(car3flg, shape = (/3*mpert*3*mpert*3*mpert/)) 4331 ddb%val(1,1:nsize,iblok) = reshape(d3cart(1,:,:,:,:,:,:), shape = (/3*mpert*3*mpert*3*mpert/)) 4332 ddb%val(2,1:nsize,iblok) = reshape(d3cart(2,:,:,:,:,:,:), shape = (/3*mpert*3*mpert*3*mpert/)) 4333 4334 ABI_FREE(d3cart) 4335 ABI_FREE(car3flg) 4336 ABI_FREE(tmpflg) 4337 ABI_FREE(tmpval) 4338 ABI_FREE(rfpert) 4339 4340 else if (ddb%typ(iblok) == BLKTYP_d3E_lw) then 4341 4342 nsize=3*mpert*3*mpert*3*mpert 4343 ABI_MALLOC(tmpflg,(3,mpert,3,mpert,3,mpert)) 4344 ABI_MALLOC(tmpval,(2,3,mpert,3,mpert,3,mpert)) 4345 4346 tmpflg(:,:,:,:,:,:) = reshape(ddb%flg(1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/)) 4347 tmpval(1,:,:,:,:,:,:) = reshape(ddb%val(1,1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/)) 4348 tmpval(2,:,:,:,:,:,:) = reshape(ddb%val(2,1:nsize,iblok), shape = (/3,mpert,3,mpert,3,mpert/)) 4349 4350 ABI_MALLOC(d3cart,(2,3,mpert,3,mpert,3,mpert)) 4351 ABI_MALLOC(car3flg,(3,mpert,3,mpert,3,mpert)) 4352 4353 call lwcart(tmpflg,car3flg,tmpval,d3cart,gprimd,mpert,natom,rprimd) 4354 4355 ddb%flg(1:nsize,iblok) = reshape(car3flg, shape = (/3*mpert*3*mpert*3*mpert/)) 4356 ddb%val(1,1:nsize,iblok) = reshape(d3cart(1,:,:,:,:,:,:), shape = (/3*mpert*3*mpert*3*mpert/)) 4357 ddb%val(2,1:nsize,iblok) = reshape(d3cart(2,:,:,:,:,:,:), shape = (/3*mpert*3*mpert*3*mpert/)) 4358 4359 ABI_FREE(d3cart) 4360 ABI_FREE(car3flg) 4361 ABI_FREE(tmpflg) 4362 ABI_FREE(tmpval) 4363 end if 4364 4365 ABI_FREE(symq) 4366 4367 end subroutine ddb_symmetrize_and_transform
m_ddb/ddb_to_dtset [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_to_dtset
FUNCTION
Initialize a dataset object from ddb. FIXME: I don't understand the goal of this routine. The dtset constructed from the DDB won't be equal to the one used to generate the DDB There's only one safe way to init dtset i.e. from file by calling the parser
INPUTS
OUTPUT
SOURCE
6009 subroutine ddb_to_dtset(comm, dtset, filename, psps) 6010 6011 !Arguments ------------------------------- 6012 integer,intent(in) :: comm 6013 type(dataset_type),intent(inout) :: dtset 6014 type(pseudopotential_type),intent(inout) :: psps 6015 ! type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 6016 character(len=*),intent(in) :: filename 6017 !Local variables ------------------------- 6018 integer :: mxnimage,unddb 6019 !integer :: ii, nn 6020 type(ddb_hdr_type) :: ddb_hdr 6021 6022 ! ************************************************************************ 6023 6024 ABI_UNUSED(psps%usepaw) 6025 6026 !Set variables 6027 mxnimage = 1 ! Only 1 image in the DDB 6028 6029 ! Must read natom from the DDB before being able to allocate some arrays needed for invars9 6030 unddb = get_unit() 6031 call ddb_hdr%open_read(filename,comm=comm) 6032 call ddb_hdr%close() 6033 !close ddb file, just want to read the headers 6034 dtset%ngfft = ddb_hdr%ngfft 6035 6036 ! Copy scalars from ddb 6037 dtset%natom = ddb_hdr%natom 6038 dtset%mband = ddb_hdr%mband 6039 dtset%nkpt = ddb_hdr%nkpt 6040 dtset%nsym = ddb_hdr%msym 6041 dtset%ntypat = ddb_hdr%ntypat 6042 dtset%nspden = ddb_hdr%nspden 6043 dtset%nspinor = ddb_hdr%nspinor 6044 dtset%nsppol = ddb_hdr%nsppol 6045 dtset%occopt = ddb_hdr%occopt 6046 dtset%usepaw = ddb_hdr%usepaw 6047 dtset%intxc = ddb_hdr%intxc 6048 dtset%ixc = ddb_hdr%ixc 6049 dtset%iscf = ddb_hdr%iscf 6050 dtset%dilatmx = ddb_hdr%dilatmx 6051 dtset%ecut = ddb_hdr%ecut 6052 dtset%ecutsm = ddb_hdr%ecutsm 6053 dtset%pawecutdg = ddb_hdr%pawecutdg 6054 dtset%kptnrm = ddb_hdr%kptnrm 6055 dtset%dfpt_sciss = ddb_hdr%dfpt_sciss 6056 dtset%tolwfr = 1.0_dp ! dummy 6057 dtset%tphysel = ddb_hdr%tphysel 6058 dtset%tsmear = ddb_hdr%tsmear 6059 6060 ! Copy arrays from ddb 6061 ABI_REMALLOC(dtset%acell_orig, (3,mxnimage)) 6062 dtset%acell_orig(1:3,1) = ddb_hdr%acell(:) 6063 6064 ABI_REMALLOC(dtset%rprim_orig, (3,3,mxnimage)) 6065 dtset%rprim_orig(1:3,1:3,1) = ddb_hdr%rprim(:,:) 6066 6067 ABI_REMALLOC(dtset%rprimd_orig, (3,3,mxnimage)) 6068 dtset%rprimd_orig(:,1,1) = ddb_hdr%rprim(:,1) * dtset%acell_orig(1,1) 6069 dtset%rprimd_orig(:,2,1) = ddb_hdr%rprim(:,2) * dtset%acell_orig(2,1) 6070 dtset%rprimd_orig(:,3,1) = ddb_hdr%rprim(:,3) * dtset%acell_orig(3,1) 6071 6072 ABI_REMALLOC(dtset%amu_orig,(dtset%ntypat,mxnimage)) 6073 dtset%amu_orig(:,1) = ddb_hdr%amu(:) 6074 6075 ABI_REMALLOC(dtset%typat, (dtset%natom)) 6076 dtset%typat(:) = ddb_hdr%typat(1:ddb_hdr%matom) 6077 6078 ABI_REMALLOC(dtset%spinat, (3,dtset%natom)) 6079 dtset%spinat(:,:) = ddb_hdr%spinat(1:3,1:ddb_hdr%matom) 6080 6081 #ifdef FC_LLVM 6082 ! LLVM 16 doesn't recognize this macro here 6083 ABI_REMALLOC(dtset%xred_orig, (3,dtset%natom,mxnimage) ) 6084 #else 6085 ABI_REMALLOC(dtset%xred_orig, (3,dtset%natom,mxnimage)) 6086 #endif 6087 dtset%xred_orig(:,:,1) = ddb_hdr%xred(1:3,1:ddb_hdr%matom) 6088 6089 ABI_REMALLOC(dtset%ziontypat, (dtset%ntypat)) 6090 dtset%ziontypat(1:ddb_hdr%mtypat) = ddb_hdr%zion(1:ddb_hdr%mtypat) 6091 6092 ABI_REMALLOC(dtset%znucl,(dtset%ntypat)) 6093 dtset%znucl(:) = ddb_hdr%znucl(1:ddb_hdr%mtypat) 6094 6095 ABI_REMALLOC(dtset%nband,(dtset%nkpt)) 6096 dtset%nband(:) = ddb_hdr%nband(1:ddb_hdr%mkpt*ddb_hdr%nsppol) 6097 6098 ABI_REMALLOC(dtset%symafm,(dtset%nsym)) 6099 dtset%symafm(:) = ddb_hdr%symafm(1:ddb_hdr%msym) 6100 6101 ABI_REMALLOC(dtset%symrel, (3,3,dtset%nsym) ) 6102 dtset%symrel(:,:,:) = ddb_hdr%symrel(1:3,1:3,1:ddb_hdr%msym) 6103 6104 ABI_REMALLOC(dtset%tnons,(3,dtset%nsym)) 6105 dtset%tnons(:,:) = ddb_hdr%tnons(1:3,1:ddb_hdr%msym) 6106 6107 ABI_REMALLOC(dtset%kpt,(3,dtset%nkpt)) 6108 dtset%kpt(:,:) = ddb_hdr%kpt(1:3,1:ddb_hdr%mkpt) 6109 6110 ABI_REMALLOC(dtset%wtk,(dtset%nkpt)) 6111 dtset%wtk(:) = ddb_hdr%wtk(1:ddb_hdr%mkpt) 6112 6113 ! GA: I had way too much problems implementing pawtab_copy. 6114 ! The script check-libpaw would report all sorts of errors. 6115 ! Therefore, I do a cheap copy here, copying only the relevant info. 6116 !call pawtab_copy(pawtab, ddb_hdr%pawtab) 6117 ! nn=size(pawtab) 6118 ! if (nn.gt.0) then 6119 ! do ii=1,nn 6120 ! pawtab(ii)%basis_size =ddb_hdr%pawtab(ii)%basis_size 6121 ! pawtab(ii)%lmn_size =ddb_hdr%pawtab(ii)%lmn_size 6122 ! pawtab(ii)%lmn2_size =ddb_hdr%pawtab(ii)%lmn2_size 6123 ! pawtab(ii)%rpaw =ddb_hdr%pawtab(ii)%rpaw 6124 ! pawtab(ii)%rshp =ddb_hdr%pawtab(ii)%rshp 6125 ! pawtab(ii)%shape_type =ddb_hdr%pawtab(ii)%shape_type 6126 ! if (allocated(pawtab(ii)%dij0)) then 6127 ! call alloc_copy(ddb_hdr%pawtab(ii)%dij0, pawtab(ii)%dij0) 6128 ! end if 6129 ! end do 6130 ! end if 6131 6132 call ddb_hdr%free() 6133 6134 end subroutine ddb_to_dtset
m_ddb/ddb_type [ Types ]
NAME
ddb_type
FUNCTION
Provides methods to extract and post-process the results in the derivative database (DDB)
SOURCE
83 type,public :: ddb_type 84 85 logical :: has_ncid_open 86 ! Is currently reading a netcdf file 87 88 integer :: iblock_d2eig_nc 89 ! Is currently reading a netcdf file 90 91 integer :: msize 92 ! Maximum size of dynamical matrices and other perturbations (ddk, dde...) 93 94 integer :: mpert 95 ! Maximum number of perturbations 96 97 integer :: nblok 98 ! Number of 2dte blocks in present object 99 100 integer :: natom 101 ! Number of atoms in the unit cell. 102 103 integer :: ntypat 104 ! Number of type of atoms. 105 106 integer :: occopt 107 ! Occupation option. 108 109 integer :: prtvol 110 ! Verbosity level. 111 112 integer :: nband 113 ! Number of bands for eigenvalues derivatives 114 ! This corresponds to d2eig arrary shape, 115 ! but the actual number of band is nband / nsppol 116 117 integer :: nkpt 118 ! Number of k-points for eigenvalues derivatives 119 120 ! GA: FIXME 121 integer :: nsppol 122 ! Number of spin components for eigenvalues derivatives 123 ! This index is absorbed into nband, to limit array ranks to 7. 124 125 integer :: current_iblok 126 ! Number of k-points for eigenvalues derivatives 127 128 ! These values are used to call the anaddb routines that don't use rprimd, gprimd. 129 real(dp) :: rprim(3,3) 130 real(dp) :: gprim(3,3) 131 real(dp) :: acell(3) 132 133 ! Many of these variables should become private so that one can refactor the ddb_t implementation 134 integer,allocatable :: flg(:,:) 135 ! flg(msize,nblok) 136 ! Flag to indicate presence of a given block 137 138 integer,allocatable :: typ(:) 139 ! typ(nblok) 140 ! Type of each block - nth-order derivatives of energy or eigenvalues. 141 ! (0 => total energy) 142 ! (1=> non-stationary block), 143 ! (2=> stationary block), 144 ! (3=> third order derivative). 145 ! (4 => first-order derivatives of total energy) 146 ! (5 => 2nd-order derivatives of eigenvalues) 147 ! (33 => long wave third order derivatives of total energy) 148 ! (85 => Molecular Berry curvature, 2nd-order derivative) 149 150 real(dp),allocatable :: amu(:) 151 ! amu(ntypat) 152 ! Mass of the atoms (atomic mass unit) 153 154 real(dp),allocatable :: qpt(:,:) 155 ! qpt(9,nblok) 156 ! q-point vector in reciprocal space (reduced lattice coordinates) for each block 157 ! Three possible phonon wavevectors can be specified for 3rd order derivatives, 158 ! but only one should be used in case of second derivative of total energy, 159 ! because we know that the second is the opposite of this value. 160 161 real(dp),allocatable :: nrm(:,:) 162 ! nrm(3,nblok) 163 ! Normalization factors of the wavevectors for each block - can be 0 to indicate a direction of approach to gamma 164 165 real(dp),allocatable :: val(:,:,:) 166 ! val(2,msize,nblok) 167 ! Values of the second energy derivatives in each block 168 169 real(dp),allocatable :: kpt(:,:) 170 ! kpt(3,nkpt) 171 ! k-point vector in reciprocal space for eigenvalues derivatives 172 173 real(dp),allocatable :: eig2dval(:,:,:,:) 174 ! eig2dval(2,msize,nband,nkpt) 175 ! Values of the second derivatives of eigenvalues 176 ! Only a single block (a single q-point) is held in memory. 177 ! Note that isppol index is wrapped into nband index. 178 179 contains 180 181 procedure :: init => ddb_init 182 ! Construct the object from the dtset. 183 184 procedure :: free => ddb_free 185 ! Free dynamic memory. 186 187 procedure :: malloc => ddb_malloc 188 ! Allocate dynamic memory 189 190 procedure :: malloc_d2eig => ddb_malloc_d2eig 191 ! Allocate dynamic memory 192 193 procedure :: copy => ddb_copy 194 ! Copy the object. 195 196 !procedure :: get_qptopt => ddb_get_qptopt 197 198 procedure :: set_qpt => ddb_set_qpt 199 ! Set the wavevector 200 201 procedure :: set_d1matr => ddb_set_d1matr 202 ! Set values for the first-order derivative matrix in tensor shape 203 204 procedure :: get_d1matr => ddb_get_d1matr 205 ! Transform the first-order derivative matrix in tensor shape 206 207 procedure :: set_d2matr => ddb_set_d2matr 208 ! Set values for the second-order derivative matrix 209 210 procedure :: get_d2matr => ddb_get_d2matr 211 ! Transform the second-order derivative matrix in tensor shape 212 213 procedure :: set_d3matr => ddb_set_d3matr 214 ! Set values for the third-order derivative matrix 215 216 procedure :: get_d3matr => ddb_get_d3matr 217 ! Transform the third-order derivative matrix in tensor shape 218 219 procedure :: get_d2eig => ddb_get_d2eig 220 ! Transform the second-order derivative matrix of eigs in tensor shape 221 222 procedure :: set_d2eig => ddb_set_d2eig 223 ! Set values for the second-order derivative matrix of eigs 224 225 procedure :: set_d2eig_reshape => ddb_set_d2eig_reshape 226 ! Set values for the second-order derivative matrix of eigs 227 ! with band index before perturbation indices 228 229 procedure :: set_gred => ddb_set_gred 230 ! Set the gradient of total energy in reduced coordinates 231 232 procedure :: set_pel => ddb_set_pel 233 ! Set the electronic polarization 234 235 procedure :: set_strten => ddb_set_strten 236 ! Set the stress tensor 237 238 procedure :: set_etotal => ddb_set_etotal 239 ! Set the total energy 240 241 procedure :: set_brav => ddb_set_brav 242 ! Set the bravais lattice. 243 244 procedure :: set_typ => ddb_set_typ 245 ! Set the typ of one block 246 247 procedure :: bcast => ddb_bcast 248 ! Broadcast the object. 249 250 procedure :: get_etotal => ddb_get_etotal 251 ! Read the GS total energy. 252 253 procedure :: get_gred => ddb_get_gred 254 ! Get the gradient of total energy in reduced coordinates 255 256 procedure :: get_pel => ddb_get_pel 257 ! Get the electronic polarization 258 259 procedure :: get_strten => ddb_get_strten 260 ! Get the stress tensor 261 262 procedure :: get_dielt_zeff => ddb_get_dielt_zeff 263 ! Reads the Dielectric Tensor and the Effective Charges 264 265 procedure :: get_dielt => ddb_get_dielt 266 ! Reads the Dielectric Tensor 267 268 procedure :: get_quadrupoles => ddb_get_quadrupoles 269 ! Reads the Quadrupoles 270 271 procedure :: get_dchidet => ddb_get_dchidet 272 ! Reads the non-linear optical susceptibility tensor and the 273 ! first-order change in the linear dielectric susceptibility 274 275 procedure :: diagoq => ddb_diagoq 276 ! Compute the phonon frequencies at the specified q-point by performing 277 ! a direct diagonalizatin of the dynamical matrix. 278 279 procedure :: get_asrq0 => ddb_get_asrq0 280 ! Return object used to enforce the acoustic sum rule 281 282 procedure :: symmetrize_and_transform => ddb_symmetrize_and_transform 283 ! Symmetrize, transform cartesian coordinates, and add missing components 284 285 procedure :: write_block_txt => ddb_write_block_txt 286 ! Writes blocks of data in the DDB in text format. 287 288 procedure :: write => ddb_write 289 ! Write the DDB file in either txt or netcdf format. 290 291 procedure :: write_txt => ddb_write_txt 292 ! Write the body of the DDB text file. 293 294 procedure :: write_nc => ddb_write_nc 295 ! Write the netcdf file (DDB.nc). 296 297 procedure :: read_block_txt => ddb_read_block_txt 298 ! Read blocks of data in the DDB. 299 300 procedure :: get_block => ddb_get_block 301 ! Finds the block containing the derivatives of the total energy. 302 303 procedure :: read_d2eig => ddb_read_d2eig 304 ! Read the next DDB block containing 2nd order derivatives of eigenvalues. 305 306 procedure :: read_d2eig_txt => ddb_read_d2eig_txt 307 ! Read the next DDB block containing 2nd order derivatives of eigenvalues. 308 309 procedure :: read_d2eig_nc => ddb_read_d2eig_nc 310 ! Read the next DDB block containing 2nd order derivatives of eigenvalues. 311 312 procedure :: write_d2eig => ddb_write_d2eig 313 ! Read the current DDB block containing 2nd order derivatives of eigenvalues. 314 315 procedure :: write_d2eig_txt => ddb_write_d2eig_txt 316 ! Read the current DDB block containing 2nd order derivatives of eigenvalues. 317 318 procedure :: write_d2eig_nc => ddb_write_d2eig_nc 319 ! Write the current DDB block containing 2nd order derivatives of eigenvalues. 320 321 procedure :: read_d0E_nc => ddb_read_d0E_nc 322 ! Read the next DDB block containing 0th order derivatives of energy. 323 324 procedure :: read_d1E_nc => ddb_read_d1E_nc 325 ! Read the next DDB block containing 1st order derivatives of energy. 326 327 procedure :: read_d2E_nc => ddb_read_d2E_nc 328 ! Read the next DDB block containing 2nd order derivatives of energy. 329 330 procedure :: read_d3E_nc => ddb_read_d3E_nc 331 ! Read the next DDB block containing 3rd order derivatives of energy. 332 333 procedure :: from_file => ddb_from_file 334 ! Construct the object from the DDB file. 335 336 procedure :: read_txt => ddb_read_txt 337 ! Construct the object from the DDB file in text format. 338 339 procedure :: read_nc => ddb_read_nc 340 ! Construct the object from the DDB file in netcdf format. 341 342 procedure :: can_merge_blocks => ddb_can_merge_blocks 343 ! Tell if two blocks can be merged 344 345 procedure :: merge_blocks => ddb_merge_blocks 346 ! Merge a block of an other ddb to the current object. 347 348 end type ddb_type 349 350 public :: ddb_to_dtset ! Transfer ddb_hdr to dtset datatype 351 public :: merge_ddb ! Read a list of ddb files and merge them into a single ddb object
m_ddb/ddb_write [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_write
FUNCTION
Write the DDB file in either txt or netcdf format.
INPUTS
ddb_hdr=ddb header object.
filename=name of the file being written (abo_DS*_DDB)
with_psps
1-> include information on pseudopoentials
0-> do not include information on pseudopoentials
comm=MPI communicator
SOURCE
4773 subroutine ddb_write(ddb, ddb_hdr, filename, with_psps, comm) 4774 4775 !Arguments ------------------------------- 4776 class(ddb_type),intent(inout) :: ddb 4777 type(ddb_hdr_type),intent(inout) :: ddb_hdr 4778 character(len=fnlen),intent(in) :: filename 4779 integer,intent(in),optional :: with_psps 4780 integer,intent(in),optional :: comm 4781 4782 !Local variables------------------------------- 4783 character(len=fnlen) :: filename_ 4784 integer :: iomode 4785 4786 ! ************************************************************************ 4787 4788 call ddb_hdr%get_iomode(filename, 2, iomode, filename_) 4789 4790 if (iomode==IO_MODE_ETSF) then 4791 call ddb%write_nc(ddb_hdr, filename_, comm=comm, with_psps=with_psps) 4792 else if (iomode==IO_MODE_FORTRAN) then 4793 call ddb%write_txt(ddb_hdr, filename_, with_psps=with_psps, comm=comm) 4794 ddb_hdr%mpert = ddb%mpert ! Text format doesnt know about mpert. 4795 end if 4796 4797 end subroutine ddb_write
m_ddb/ddb_write_block_txt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_write_block_txt
FUNCTION
This routine writes blocks of data in the DDB in text format.
INPUTS
choice= (2 => write), (3 => write minimal info ) mpert =maximum number of ipert msize=maximum size of the arrays flags and values nunit=unit number for the data block file
OUTPUT
(see side effects)
SIDE EFFECTS
Input/Output ddb = ddb block datastructure ddb%typ=type of the block: 0 => total energy 1 => second-order energy derivatives, non-stationary block 2 => second-order energy derivatives, stationary block 3 => third-order energy derivatives 4 => first-order energy derivatives: forces, stresses and polarization 5 => second-order eigenvalue derivatives ddb%flg(msize)=flag for every matrix element (0=> the element is not in the data block), (1=> the element is in the data blok) ddb%qpt(9)=wavevector of the perturbation(s). The elements from 1 to 3 are used if we are dealing with the 2nd derivative of total energy (only one wavevector), while all elements are used in case of a third order derivative of total energy (three wavevector could be present) ddb%nrm(3)=normalization factors for the three allowed wavevectors. ddb%val(2,msize)=real(dp), complex, value of the matrix elements that are present in the data block blkval2(2,msize,mband,nkpt) = value of the matrix elements that are present in a block of EIGR2D/EIGI2D
NOTES
only executed by one processor.
SOURCE
4587 subroutine ddb_write_block_txt(ddb,iblok,choice,mband,mpert,msize,nkpt,nunit,& 4588 blkval2,kpt) !optional 4589 4590 !Arguments ------------------------------- 4591 !scalars 4592 integer,intent(in) :: choice,mband,mpert,msize,nkpt,nunit 4593 integer,intent(in) :: iblok 4594 class(ddb_type),intent(in) :: ddb 4595 !arrays 4596 real(dp),intent(in),optional :: kpt(3,nkpt) 4597 real(dp),intent(in),optional :: blkval2(2,msize,mband,nkpt) 4598 4599 !Local variables ------------------------- 4600 !scalars 4601 integer :: iband,idir1,idir2,idir3,ii,ikpt,ipert1,ipert2,ipert3 4602 integer :: nelmts 4603 logical :: eig2d_ 4604 4605 ! ********************************************************************* 4606 4607 ! GA: Remove choice option. 4608 ! choice=3 is used to write summary info at the end of the DDB. 4609 ! With MG and MV, we agreed that it could be removed. 4610 4611 ! GA: Remove arguments: mband, blkval2, kpt 4612 ! This feature of writing eigenvalues 2nd deriv is no longer used 4613 ! (see 80_tdep/m_tdep_abitypes.F90) 4614 ! With FB, we agreed that it could be removed. 4615 4616 eig2d_ = .false. 4617 if(present(blkval2).and.present(kpt)) eig2d_ = .true. 4618 4619 4620 ! Count the number of elements 4621 nelmts=0 4622 do ii=1,msize 4623 if(ddb%flg(ii,iblok)==1)nelmts=nelmts+1 4624 end do 4625 4626 ! Write the block type and number of elements 4627 write(nunit,*)' ' 4628 if (ddb%typ(iblok) == BLKTYP_d0E_xx) then 4629 write(nunit, '(a,i12)' )' Total energy - # elements :',nelmts 4630 else if (ddb%typ(iblok)==BLKTYP_d2E_ns) then 4631 write(nunit, '(a,i12)' )' 2nd derivatives (non-stat.) - # elements :',nelmts 4632 else if(ddb%typ(iblok)==BLKTYP_d2E_st) then 4633 write(nunit, '(a,i12)' )' 2nd derivatives (stationary) - # elements :',nelmts 4634 else if(ddb%typ(iblok)==BLKTYP_d2E_mbc) then 4635 write(nunit, '(a,i12)' )' 2nd derivatives (MBC) - # elements :',nelmts 4636 else if(ddb%typ(iblok)==BLKTYP_d3E_xx) then 4637 write(nunit, '(a,i12)' )' 3rd derivatives - # elements :',nelmts 4638 else if (ddb%typ(iblok) == BLKTYP_d1E_xx) then 4639 write(nunit, '(a,i12)' )' 1st derivatives - # elements :',nelmts 4640 else if (ddb%typ(iblok) == BLKTYP_d2eig_re) then 4641 write(nunit, '(a,i12)' )' 2nd eigenvalue derivatives - # elements :',nelmts 4642 else if(ddb%typ(iblok)==BLKTYP_d3E_lw) then 4643 write(nunit, '(a,i12)' )' 3rd derivatives (long wave) - # elements :',nelmts 4644 end if 4645 4646 ! Write the 2nd derivative block 4647 if (is_type_d2E(ddb%typ(iblok))) then 4648 4649 ! Write the phonon wavevector 4650 write(nunit, '(a,3es16.8,f6.1)' )' qpt',(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok) 4651 4652 ! Write the matrix elements 4653 if(choice==2)then 4654 ii=0 4655 do ipert2=1,mpert 4656 do idir2=1,3 4657 do ipert1=1,mpert 4658 do idir1=1,3 4659 ii=ii+1 4660 if(ddb%flg(ii,iblok)==1)then 4661 write(nunit,'(4i4,2d22.14)') idir1, ipert1, idir2, ipert2, & 4662 ddb%val(1,ii,iblok),ddb%val(2,ii,iblok) 4663 end if 4664 end do 4665 end do 4666 end do 4667 end do 4668 end if 4669 4670 4671 else if (is_type_d3E(ddb%typ(iblok))) then 4672 ! Write the 3rd derivative block 4673 4674 ! Write the phonon wavevectors 4675 write(nunit, '(a,3es16.8,f6.1)' )' qpt',(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok) 4676 write(nunit, '(a,3es16.8,f6.1)' )' ',(ddb%qpt(ii,iblok),ii=4,6),ddb%nrm(2,iblok) 4677 write(nunit, '(a,3es16.8,f6.1)' )' ',(ddb%qpt(ii,iblok),ii=7,9),ddb%nrm(3,iblok) 4678 4679 ! Write the matrix elements 4680 if(choice==2)then 4681 ii=0 4682 do ipert3=1,mpert 4683 do idir3=1,3 4684 do ipert2=1,mpert 4685 do idir2=1,3 4686 do ipert1=1,mpert 4687 do idir1=1,3 4688 ii=ii+1 4689 if(ddb%flg(ii,iblok)==1)then 4690 write(nunit, '(6i6,2d22.14)' )& 4691 idir1,ipert1,idir2,ipert2,idir3,ipert3,ddb%val(1,ii,iblok),ddb%val(2,ii,iblok) 4692 end if 4693 end do 4694 end do 4695 end do 4696 end do 4697 end do 4698 end do 4699 end if 4700 4701 4702 else if (is_type_d0E(ddb%typ(iblok))) then 4703 ! Write total energy 4704 if (choice == 2) write(nunit,'(2d22.14)')ddb%val(1,1,iblok),ddb%val(2,1,iblok) 4705 4706 else if (is_type_d1E(ddb%typ(iblok))) then 4707 ! Write the 1st derivative blok 4708 if (choice == 2) then 4709 ii = 0 4710 do ipert1 = 1, mpert 4711 do idir1 = 1, 3 4712 ii = ii + 1 4713 if (ddb%flg(ii,iblok) == 1) then 4714 write(nunit,'(2i6,2d22.14)')idir1,ipert1,ddb%val(1,ii,iblok),ddb%val(2,ii,iblok) 4715 end if 4716 end do 4717 end do 4718 end if 4719 4720 else if (is_type_d2eig(ddb%typ(iblok))) then 4721 ! Write the phonon wavevector 4722 write(nunit, '(a,3es16.8,f6.1)' )' qpt',(ddb%qpt(ii,iblok),ii=1,3),ddb%nrm(1,iblok) 4723 ! Write the matrix elements 4724 ! GA: Note that isppol is invisible here. It is simply marked as more bands. 4725 ! To be changed in a future version of text format. 4726 if(choice==2)then 4727 if (eig2d_) then 4728 do ikpt=1,nkpt 4729 write(nunit,'(a,3es16.8)')' K-point:',(kpt(ii,ikpt),ii=1,3) 4730 do iband=1,mband 4731 write(nunit,'(a,i3)')' Band:',iband 4732 ii=0 4733 do ipert2=1,mpert 4734 do idir2=1,3 4735 do ipert1=1,mpert 4736 do idir1=1,3 4737 ii=ii+1 4738 if(ddb%flg(ii,iblok)==1)then 4739 write(nunit,'(4i4,2d22.14)')idir1,ipert1,idir2,ipert2,blkval2(1,ii,iband,ikpt),blkval2(2,ii,iband,ikpt) 4740 end if 4741 end do !idir1 4742 end do !ipert1 4743 end do !idir2 4744 end do !ipert2 4745 end do !iband 4746 end do !ikpt 4747 end if !eig2d_ 4748 end if !choice 4749 end if !ddb%typ(iblok) 4750 4751 end subroutine ddb_write_block_txt
m_ddb/ddb_write_d2eig [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_write_d2eig
FUNCTION
Write the current eig2d data as the next block in the ddb file.
INPUTS
ddb_hdr=ddb header object. unddb=unit of the open ddb file in text format or netcdf identifier.
SOURCE
4865 subroutine ddb_write_d2eig(ddb, ddb_hdr, iblok, comm) 4866 !Arguments ------------------------------- 4867 class(ddb_type),intent(inout) :: ddb 4868 type(ddb_hdr_type),intent(inout) :: ddb_hdr 4869 integer,intent(in) :: iblok 4870 integer,intent(in),optional :: comm 4871 4872 !Local variables ------------------------- 4873 !scalars 4874 integer,parameter :: master=0 4875 character(len=500) :: msg 4876 4877 ! ************************************************************************ 4878 4879 if (present(comm)) then 4880 if (xmpi_comm_rank(comm) /= master) return 4881 end if 4882 4883 if (ddb_hdr%has_open_file_nc) then 4884 4885 call ddb%write_d2eig_nc(ddb_hdr%ncid, iblok) 4886 4887 else if (ddb_hdr%has_open_file_txt) then 4888 4889 call ddb%write_d2eig_txt(ddb_hdr%unddb, iblok) 4890 4891 else 4892 write(msg, '(3a)' )& 4893 ! File has not been opened by ddb_hdr 4894 'Attempting to write into unopen DDB file.',ch10,& 4895 'Action: contact Abinit group.' 4896 ABI_ERROR(msg) 4897 end if 4898 4899 end subroutine ddb_write_d2eig
m_ddb/ddb_write_d2eig_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_write_d2eig_nc
FUNCTION
Write the current d2eig data in the ddb netcdf file.
INPUTS
iblok=index of the eig2d block within the d2eig subgroup. ncid=netcdf identifier of a file open in writing mode. comm=MPI communicator.
SOURCE
4918 subroutine ddb_write_d2eig_nc(ddb, ncid, iblok, comm) 4919 !Arguments ------------------------------- 4920 class(ddb_type),intent(inout) :: ddb 4921 integer,intent(in) :: ncid 4922 integer,intent(in) :: iblok 4923 integer,intent(in),optional :: comm 4924 4925 !Local variables ------------------------- 4926 !scalars 4927 integer,parameter :: master=0 4928 integer :: iband, jband, bandshift, isppol, mband 4929 integer :: ikpt, ipert1, idir1, ipert2, idir2, ii 4930 integer :: ncid_d2eig, ncerr 4931 real(dp) :: qpt(3) 4932 real(dp), allocatable :: matrix_d2eig(:,:,:,:,:,:,:) 4933 real(dp), allocatable :: matrix_d2eig_isppol(:,:,:,:,:,:,:) 4934 integer, allocatable :: flg_d2eig(:,:,:,:) 4935 4936 ! ************************************************************************ 4937 4938 4939 if (present(comm)) then 4940 if (xmpi_comm_rank(comm) /= master) return 4941 end if 4942 4943 ncid_d2eig = nctk_idgroup(ncid, 'd2eig') 4944 4945 qpt(1:3) = ddb%qpt(1:3,iblok) 4946 ncerr = nf90_put_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 4947 'reduced_coordinates_of_qpoints'),& 4948 qpt,& 4949 start=[1,iblok]) 4950 NCF_CHECK(ncerr) 4951 ncerr = nf90_put_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 4952 'qpoints_normalization'),& 4953 ddb%nrm(1,iblok),& 4954 start=[iblok]) 4955 NCF_CHECK(ncerr) 4956 4957 ! GA: Here we assume that all blocks are d2eig blocks. 4958 ! Otherwise, iblok is only the 'local' iblok index 4959 ! and we would need to figure out the corresponding 'global' index 4960 call ddb%get_d2eig(matrix_d2eig, flg_d2eig, iblok) 4961 4962 mband = ddb%nband / ddb%nsppol 4963 ABI_MALLOC(matrix_d2eig_isppol, (2,3,ddb%mpert,3,ddb%mpert,mband,ddb%nkpt)) 4964 4965 ! Loop over spin index 4966 do isppol=1,ddb%nsppol 4967 4968 bandshift = (isppol - 1) * mband 4969 4970 do iband=1,mband 4971 jband = bandshift + iband 4972 4973 do ikpt=1,ddb%nkpt 4974 do ipert1=1,ddb%mpert 4975 do idir1=1,3 4976 do ipert2=1,ddb%mpert 4977 do idir2=1,3 4978 do ii=1,2 4979 matrix_d2eig_isppol(ii,idir2,ipert2,idir1,ipert1,iband,ikpt)=& 4980 matrix_d2eig(ii,idir2,ipert2,idir1,ipert1,jband,ikpt) 4981 end do 4982 end do 4983 end do 4984 end do 4985 end do 4986 end do 4987 4988 4989 end do ! iband 4990 4991 ncerr = nf90_put_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 4992 'matrix_values'),& 4993 matrix_d2eig_isppol,& 4994 start=[1,1,1,1,1,1,1,isppol,iblok]) 4995 !count=[2,3,ddb%mpert,3,ddb%mpert,ddb%nband,ddb%nkpt,1,1]) 4996 !count=[2,3,ddb%mpert,3,ddb%mpert,ddb%nkpt,ddb%nband,1,1]) 4997 NCF_CHECK(ncerr) 4998 4999 end do ! isppol 5000 5001 ncerr = nf90_put_var(ncid_d2eig, nctk_idname(ncid_d2eig,& 5002 'matrix_mask'),& 5003 flg_d2eig,& 5004 start=[1,1,1,1,iblok]) 5005 NCF_CHECK(ncerr) 5006 5007 ABI_SFREE(matrix_d2eig_isppol) 5008 ABI_SFREE(matrix_d2eig) 5009 ABI_SFREE(flg_d2eig) 5010 5011 end subroutine ddb_write_d2eig_nc
m_ddb/ddb_write_d2eig_txt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_write_d2eig_txt
FUNCTION
Write the eig2d data as the next block in text file format.
INPUTS
ddb_hdr=ddb header object. unddb=unit of the open ddb file in text format.
SOURCE
5031 subroutine ddb_write_d2eig_txt(ddb, unddb, iblok) 5032 !Arguments ------------------------------- 5033 class(ddb_type),intent(in) :: ddb 5034 integer,intent(in) :: unddb 5035 integer,intent(in) :: iblok 5036 5037 !Local variables ------------------------- 5038 !scalars 5039 integer,parameter :: iblok_eig2d=1 5040 integer,parameter :: choice=2 5041 5042 ! ************************************************************************ 5043 5044 ! GA: This routine is redundant with outbsd. 5045 ! The present implementation should replace outbsd. 5046 5047 call ddb%write_block_txt(iblok,choice,ddb%nband,ddb%mpert,ddb%msize,ddb%nkpt,unddb,& 5048 ddb%eig2dval(:,:,:,:), ddb%kpt(:,:)) 5049 5050 end subroutine ddb_write_d2eig_txt
m_ddb/ddb_write_nc [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_write_nc
FUNCTION
Write the netndf file DDB.nc in the format of version 20230219.
INPUTS
ddb_hdr=ddb header object with open file. filename=DDB filename. comm=MPI communicator. with_psps 1-> include information on pseudopoentials 0-> do not include information on pseudopoentials
SOURCE
5072 subroutine ddb_write_nc(ddb, ddb_hdr, filename, comm, with_psps) 5073 5074 !Arguments ------------------------------- 5075 class(ddb_type),intent(inout) :: ddb 5076 type(ddb_hdr_type),intent(inout) :: ddb_hdr 5077 character(len=*),intent(in) :: filename 5078 integer,intent(in),optional :: comm 5079 integer,intent(in),optional :: with_psps 5080 5081 !Local variables ------------------------- 5082 !scalars 5083 integer,parameter :: master=0 5084 integer :: ncid, ncerr, ncid_d0E, ncid_d1E, ncid_d2E, ncid_d3E, ncid_d2eig 5085 integer :: ii,iblok,iblok_d0E,iblok_d1E,iblok_d2E,iblok_d3E,iblok_d2eig 5086 !arrays 5087 integer,allocatable :: flg_d1E(:,:) 5088 integer,allocatable :: flg_d2E(:,:,:,:) 5089 integer,allocatable :: flg_d3E(:,:,:,:,:,:) 5090 real(dp) :: qpt(3), qpts(3,3), nrms(3) 5091 real(dp),allocatable :: matrix_d1E(:,:,:) 5092 real(dp),allocatable :: matrix_d2E(:,:,:,:,:) 5093 real(dp),allocatable :: matrix_d3E(:,:,:,:,:,:,:) 5094 5095 ! ************************************************************************ 5096 5097 if (present(comm)) then 5098 if (xmpi_comm_rank(comm) /= master) return 5099 end if 5100 5101 #ifdef HAVE_NETCDF 5102 5103 ! ===================== 5104 ! Header and dimensions 5105 ! ===================== 5106 ! Copy types and dimensions into the header 5107 ddb_hdr%mpert = ddb%mpert 5108 call ddb_hdr%set_typ(ddb%nblok, ddb%typ) 5109 5110 call ddb_hdr%open_write_nc(filename, with_psps=with_psps) 5111 ncid = ddb_hdr%ncid 5112 5113 ! Get all group id 5114 ncid_d0E = nctk_idgroup(ncid, 'd0E') 5115 ncid_d1E = nctk_idgroup(ncid, 'd1E') 5116 ncid_d2E = nctk_idgroup(ncid, 'd2E') 5117 ncid_d3E = nctk_idgroup(ncid, 'd3E') 5118 ncid_d2eig = nctk_idgroup(ncid, 'd2eig') 5119 5120 ! ============================= 5121 ! Loop over block to be written 5122 ! ============================= 5123 5124 iblok_d0E = 0; iblok_d1E = 0; iblok_d2E = 0; iblok_d3E = 0; iblok_d2eig = 0 5125 5126 do iblok=1,ddb%nblok 5127 5128 ! ------------------------ 5129 ! Zeroth-order derivatives 5130 ! ------------------------ 5131 if (is_type_d0E(ddb%typ(iblok))) then 5132 5133 iblok_d0E = iblok_d0E + 1 5134 5135 ncerr = nf90_put_var(ncid_d0E, nctk_idname(ncid_d0E,& 5136 'matrix_values'),& 5137 ddb%val(1,1,iblok),& 5138 start=[iblok_d0E]) 5139 NCF_CHECK(ncerr) 5140 5141 ncerr = nf90_put_var(ncid_d0E, nctk_idname(ncid_d0E,& 5142 'matrix_mask'),& 5143 ddb%flg(1,iblok),& 5144 start=[iblok_d0E]) 5145 NCF_CHECK(ncerr) 5146 5147 ! ----------------------- 5148 ! First-order derivatives 5149 ! ----------------------- 5150 else if (is_type_d1E(ddb%typ(iblok))) then 5151 5152 iblok_d1E = iblok_d1E + 1 5153 5154 call ddb%get_d1matr(iblok, matrix_d1E, flg_d1E) 5155 5156 ncerr = nf90_put_var(ncid_d1E, nctk_idname(ncid_d1E,& 5157 'matrix_values'),& 5158 matrix_d1E,& 5159 start=[1,1,1,iblok_d1E]) 5160 NCF_CHECK(ncerr) 5161 5162 ncerr = nf90_put_var(ncid_d1E, nctk_idname(ncid_d1E,& 5163 'matrix_mask'),& 5164 flg_d1E,& 5165 start=[1,1,iblok_d1E]) 5166 NCF_CHECK(ncerr) 5167 ABI_SFREE(matrix_d1E) 5168 ABI_SFREE(flg_d1E) 5169 5170 ! ------------------------ 5171 ! Second-order derivatives 5172 ! ------------------------ 5173 else if (is_type_d2E(ddb%typ(iblok))) then 5174 5175 iblok_d2E = iblok_d2E + 1 5176 5177 do ii=1,3 5178 qpt(ii) = ddb%qpt(ii,iblok) 5179 end do 5180 ncerr = nf90_put_var(ncid_d2E, nctk_idname(ncid_d2E,& 5181 'reduced_coordinates_of_qpoints'),& 5182 qpt,& 5183 start=[1,iblok_d2E]) 5184 NCF_CHECK(ncerr) 5185 5186 ncerr = nf90_put_var(ncid_d2E, nctk_idname(ncid_d2E,& 5187 'qpoints_normalization'),& 5188 (ddb%nrm(1,iblok)),& 5189 start=[iblok_d2E]) 5190 NCF_CHECK(ncerr) 5191 5192 call ddb%get_d2matr(iblok, matrix_d2E, flg_d2E) 5193 5194 ncerr = nf90_put_var(ncid_d2E, nctk_idname(ncid_d2E,& 5195 'matrix_values'),& 5196 matrix_d2E,& 5197 start=[1,1,1,1,1,iblok_d2E]) 5198 NCF_CHECK(ncerr) 5199 5200 ncerr = nf90_put_var(ncid_d2E, nctk_idname(ncid_d2E,& 5201 'matrix_mask'),& 5202 flg_d2E,& 5203 start=[1,1,1,1,iblok_d2E]) 5204 NCF_CHECK(ncerr) 5205 ABI_SFREE(matrix_d2E) 5206 ABI_SFREE(flg_d2E) 5207 5208 ! ----------------------- 5209 ! Third-order derivatives 5210 ! ----------------------- 5211 else if (is_type_d3E(ddb%typ(iblok))) then 5212 5213 iblok_d3E = iblok_d3E + 1 5214 5215 do ii=1,3 5216 nrms(ii) = ddb%nrm(ii,iblok) 5217 qpts(1,ii) = ddb%qpt(ii,iblok) 5218 qpts(2,ii) = ddb%qpt(ii+3,iblok) 5219 qpts(3,ii) = ddb%qpt(ii+6,iblok) 5220 end do 5221 5222 ncerr = nf90_put_var(ncid_d3E, nctk_idname(ncid_d3E,& 5223 'reduced_coordinates_of_qpoints'),& 5224 qpts,& 5225 start=[1,1,iblok_d3E]) 5226 NCF_CHECK(ncerr) 5227 5228 ncerr = nf90_put_var(ncid_d3E, nctk_idname(ncid_d3E,& 5229 'qpoints_normalization'),& 5230 nrms,& 5231 start=[1,iblok_d3E]) 5232 NCF_CHECK(ncerr) 5233 5234 call ddb%get_d3matr(iblok, matrix_d3E, flg_d3E) 5235 5236 ncerr = nf90_put_var(ncid_d3E, nctk_idname(ncid_d3E,& 5237 'matrix_values'),& 5238 matrix_d3E,& 5239 start=[1,1,1,1,1,1,1,iblok_d3E]) 5240 NCF_CHECK(ncerr) 5241 5242 ncerr = nf90_put_var(ncid_d3E, nctk_idname(ncid_d3E,& 5243 'matrix_mask'),& 5244 flg_d3E,& 5245 start=[1,1,1,1,1,1,iblok_d3E]) 5246 NCF_CHECK(ncerr) 5247 5248 ABI_SFREE(matrix_d3E) 5249 ABI_SFREE(flg_d3E) 5250 5251 ! --------------------------------------- 5252 ! Second-order derivatives of eigenvalues 5253 ! --------------------------------------- 5254 else if (is_type_d2eig(ddb%typ(iblok))) then 5255 5256 iblok_d2eig = iblok_d2eig + 1 5257 5258 call ddb%write_d2eig_nc(ncid_d2eig, iblok_d2eig) 5259 5260 end if 5261 end do 5262 5263 #else 5264 ABI_ERROR("NETCDF support required to write DDB.nc file.") 5265 #endif 5266 5267 end subroutine ddb_write_nc
m_ddb/ddb_write_txt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
ddb_write_txt
FUNCTION
Write the DDB file in text format.
INPUTS
ddb_hdr=ddb header object.
filename=name of the file being written (abo_DS*_DDB)
with_psps
1-> include information on pseudopoentials
0-> do not include information on pseudopoentials
SOURCE
4818 subroutine ddb_write_txt(ddb, ddb_hdr, filename, with_psps, comm) 4819 4820 !Arguments ------------------------------- 4821 class(ddb_type),intent(inout) :: ddb 4822 type(ddb_hdr_type),intent(inout) :: ddb_hdr 4823 character(len=*),intent(in) :: filename 4824 integer,intent(in),optional :: with_psps 4825 integer,intent(in),optional :: comm 4826 4827 !Local variables ------------------------- 4828 !scalars 4829 integer :: iblok 4830 integer,parameter :: master=0, choice=2 4831 4832 ! ************************************************************************ 4833 4834 if (present(comm)) then 4835 if (xmpi_comm_rank(comm) /= master) return 4836 end if 4837 4838 call ddb_hdr%open_write_txt(filename, with_psps) 4839 4840 do iblok=1,ddb%nblok 4841 call ddb%write_block_txt(iblok,choice,1,ddb%mpert,ddb%msize,ddb_hdr%nkpt,ddb_hdr%unddb) 4842 end do 4843 4844 call ddb_hdr%close() 4845 4846 end subroutine ddb_write_txt
m_ddb/dtchi [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
dtchi
FUNCTION
Reads the non-linear optical susceptibility tensor and the first-order change in the linear dielectric susceptibility induced by an atomic displacement in the Gamma Block coming from the Derivative Data Base (third-order derivatives).
INPUTS
blkval(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies
natom= number of atoms in unit cell
mpert =maximum number of ipert
ramansr= if /= 0, impose sum rule on first-order derivatives
of the electronic susceptibility with respect
to atomic displacements
nlflag= if =3, only the non-linear optical susceptibilities is computed
OUTPUT
dchide(3,3,3) = non-linear optical coefficients dchidt(natom,3,3,3) = first-order change of the electronic dielectric tensor induced by an individual atomic displacement
SOURCE
3333 subroutine dtchi(blkval,dchide,dchidt,mpert,natom,ramansr,nlflag) 3334 3335 !Arguments ------------------------------- 3336 !scalars 3337 integer,intent(in) :: mpert,natom,ramansr,nlflag 3338 !arrays 3339 real(dp),intent(in) :: blkval(2,3*mpert*3*mpert*3*mpert) 3340 real(dp),intent(out) :: dchide(3,3,3),dchidt(natom,3,3,3) 3341 3342 !Local variables ------------------------- 3343 !scalars 3344 integer :: depl,elfd1,elfd2,elfd3,iatom,ivoigt 3345 logical :: iwrite 3346 real(dp) :: wttot 3347 !arrays 3348 integer :: voigtindex(6,2) 3349 real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert),dvoigt(3,6),sumrule(3,3,3) 3350 real(dp) :: wghtat(natom) 3351 3352 ! ********************************************************************* 3353 3354 d3cart(1,:,:,:,:,:,:) = reshape(blkval(1,:),shape = (/3,mpert,3,mpert,3,mpert/)) 3355 d3cart(2,:,:,:,:,:,:) = reshape(blkval(2,:),shape = (/3,mpert,3,mpert,3,mpert/)) 3356 3357 ! Extraction of non-linear optical coefficients 3358 do elfd1 = 1,3 3359 do elfd2 = 1,3 3360 do elfd3 = 1,3 3361 dchide(elfd1,elfd2,elfd3) = d3cart(1,elfd1,natom+2,elfd2,natom+2,elfd3,natom+2) 3362 end do 3363 end do 3364 end do 3365 3366 ! Transform to Voigt notations 3367 voigtindex(:,1) = (/1,2,3,2,1,1/) 3368 voigtindex(:,2) = (/1,2,3,3,3,2/) 3369 do ivoigt = 1, 6 3370 elfd2 = voigtindex(ivoigt,1) 3371 elfd3 = voigtindex(ivoigt,2) 3372 do elfd1 = 1, 3 3373 dvoigt(elfd1,ivoigt) = 0.5_dp*(dchide(elfd1,elfd2,elfd3) + dchide(elfd1,elfd3,elfd2)) 3374 end do 3375 end do 3376 3377 ! Transform to pm/V 3378 dvoigt(:,:) = dvoigt(:,:)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 3379 3380 ! Extraction of $\frac{d \chi}{d \tau}$ 3381 if (nlflag < 3) then 3382 do iatom = 1, natom 3383 do depl = 1,3 3384 do elfd1 = 1,3 3385 do elfd2 = 1,3 3386 dchidt(iatom,depl,elfd1,elfd2) = d3cart(1,depl,iatom,elfd1,natom+2,elfd2,natom+2) 3387 end do 3388 end do 3389 end do 3390 end do 3391 end if 3392 3393 wghtat(:) = zero 3394 if (ramansr == 1) then 3395 wghtat(:) = one/dble(natom) 3396 3397 else if (ramansr == 2) then 3398 3399 wttot = zero 3400 do iatom = 1, natom 3401 do depl = 1,3 3402 do elfd1 = 1,3 3403 do elfd2 = 1,3 3404 wghtat(iatom) = wghtat(iatom) + abs(dchidt(iatom,depl,elfd1,elfd2)) 3405 end do 3406 end do 3407 end do 3408 wttot = wttot + wghtat(iatom) 3409 end do 3410 3411 wghtat(:) = wghtat(:)/wttot 3412 end if 3413 3414 iwrite = ab_out > 0 3415 3416 if (iwrite) then 3417 write(ab_out,*)ch10 3418 write(ab_out,*)'Non-linear optical coefficients d (pm/V)' 3419 write(ab_out,'(6f12.6)')dvoigt(1,:) 3420 write(ab_out,'(6f12.6)')dvoigt(2,:) 3421 write(ab_out,'(6f12.6)')dvoigt(3,:) 3422 end if 3423 3424 if (ramansr /= 0) then 3425 if (iwrite) then 3426 write(ab_out,*)ch10 3427 write(ab_out,*)'The violation of the Raman sum rule' 3428 write(ab_out,*)'by the first-order electronic dielectric tensors ','is as follows' 3429 write(ab_out,*)' atom' 3430 write(ab_out,*)' displacement' 3431 end if 3432 3433 sumrule(:,:,:) = zero 3434 do elfd2 = 1,3 3435 do elfd1 = 1,3 3436 do depl = 1,3 3437 do iatom = 1, natom 3438 sumrule(depl,elfd1,elfd2) = sumrule(depl,elfd1,elfd2) + dchidt(iatom,depl,elfd1,elfd2) 3439 end do 3440 do iatom = 1, natom 3441 dchidt(iatom,depl,elfd1,elfd2) = dchidt(iatom,depl,elfd1,elfd2) - wghtat(iatom)*sumrule(depl,elfd1,elfd2) 3442 end do 3443 end do 3444 end do 3445 end do 3446 3447 if (iwrite) then 3448 do depl = 1,3 3449 write(ab_out,'(6x,i2,3(3x,f16.9))') depl,sumrule(depl,1,1:3) 3450 write(ab_out,'(8x,3(3x,f16.9))') sumrule(depl,2,1:3) 3451 write(ab_out,'(8x,3(3x,f16.9))') sumrule(depl,3,1:3) 3452 write(ab_out,*) 3453 end do 3454 end if 3455 end if ! ramansr 3456 3457 if (nlflag < 3) then 3458 if (iwrite) then 3459 write(ab_out,*)ch10 3460 write(ab_out,*)' First-order change in the electronic dielectric ' 3461 write(ab_out,*)' susceptibility tensor (Bohr^-1)' 3462 write(ab_out,*)' induced by an atomic displacement' 3463 if (ramansr /= 0) then 3464 write(ab_out,*)' (after imposing the sum over all atoms to vanish)' 3465 end if 3466 write(ab_out,*)' atom displacement' 3467 3468 do iatom = 1,natom 3469 do depl = 1,3 3470 write(ab_out,'(1x,i4,9x,i2,3(3x,f16.9))')iatom,depl,dchidt(iatom,depl,1,:) 3471 write(ab_out,'(16x,3(3x,f16.9))')dchidt(iatom,depl,2,:) 3472 write(ab_out,'(16x,3(3x,f16.9))')dchidt(iatom,depl,3,:) 3473 end do 3474 3475 write(ab_out,*) 3476 end do 3477 end if 3478 end if 3479 3480 end subroutine dtchi
m_ddb/dtech9 [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
dtech9
FUNCTION
Reads the Dielectric Tensor and the Effective Charges in the Gamma Block coming from the Derivative Data Base.
INPUTS
natom= number of atoms in unit cell iblok= index of the Gamma block mpert =maximum number of ipert nblok= number of blocks in the DDB blkval(2,3*mpert*3*mpert,nblok)= dynamical matrices In our case, the nblok is restricted to iblok [unit]=Output unit number
OUTPUT
zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement. Note the following convention: zeff(electric field direction, atomic direction, atom index) dielt(3,3)=dielectric tensor
SOURCE
3250 subroutine dtech9(blkval,dielt,iblok,mpert,natom,nblok,zeff,unit) 3251 3252 !Arguments ------------------------------- 3253 !scalars 3254 integer,intent(in) :: iblok,mpert,natom,nblok 3255 integer,intent(in),optional :: unit 3256 !arrays 3257 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok) 3258 real(dp),intent(out) :: dielt(3,3),zeff(3,3,natom) 3259 3260 !Local variables ------------------------- 3261 !scalars 3262 integer :: depl,elec,elec1,elec2,iatom, unt 3263 character(len=1000) :: msg 3264 3265 ! ********************************************************************* 3266 3267 unt = std_out; if (present(unit)) unt = unit 3268 3269 ! Extraction of effectives charges 3270 do iatom=1,natom 3271 do elec=1,3 3272 do depl=1,3 3273 zeff(elec,depl,iatom)=0.5*& 3274 (blkval(1,depl,iatom,elec,natom+2,iblok)+& 3275 blkval(1,elec,natom+2,depl,iatom,iblok)) 3276 end do 3277 end do 3278 end do 3279 3280 ! Extraction of dielectric tensor 3281 do elec1=1,3 3282 do elec2=1,3 3283 dielt(elec1,elec2)=blkval(1,elec1,natom+2,elec2,natom+2,iblok) 3284 end do 3285 end do 3286 3287 write(msg,'(a,3es16.6,3es16.6,3es16.6)' )' Dielectric Tensor ',& 3288 dielt(1,1),dielt(1,2),dielt(1,3),& 3289 dielt(2,1),dielt(2,2),dielt(2,3),& 3290 dielt(3,1),dielt(3,2),dielt(3,3) 3291 call wrtout(unt, msg) 3292 3293 call wrtout(unt, ' Effectives Charges ') 3294 do iatom=1,natom 3295 write(msg,'(a,i4,3es16.6,3es16.6,3es16.6)' )' atom ',iatom,& 3296 zeff(1,1,iatom),zeff(1,2,iatom),zeff(1,3,iatom),& 3297 zeff(2,1,iatom),zeff(2,2,iatom),zeff(2,3,iatom),& 3298 zeff(3,1,iatom),zeff(3,2,iatom),zeff(3,3,iatom) 3299 call wrtout(unt, msg) 3300 end do 3301 3302 end subroutine dtech9
m_ddb/dtqdrp [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
dtqdrp
FUNCTION
Reads the Dynamical Quadrupole or the P^(1) Tensor in the Gamma Block coming from the Derivative Data Base (long wave third-order derivatives).
INPUTS
blkval(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies
ddb_version = 8 digit integer giving date. To mantain compatibility with olderDDB files.
lwsym = 0 do not symmetrize the tensor wrt efield and qvec derivative
|-> 1st gradient of polarization response to atomic displacement
= 1 symmetrize the tensor wrt efield and qvec derivative
|-> dynamic quadrupoles
natom= number of atoms in unit cell
mpert =maximum number of ipert
OUTPUT
lwtens(3,3,3,natom) = Dynamical Quadrupoles or P^(1) tensor
SOURCE
6582 subroutine dtqdrp(blkval,ddb_version,lwsym,mpert,natom,lwtens) 6583 6584 !Arguments ------------------------------- 6585 !scalars 6586 integer,intent(in) :: ddb_version,lwsym,mpert,natom 6587 !arrays 6588 real(dp),intent(in) :: blkval(2,3*mpert*3*mpert*3*mpert) 6589 real(dp),intent(out) :: lwtens(3,3,3,natom) 6590 6591 !Local variables ------------------------- 6592 !scalars 6593 integer,parameter :: cvrsio8=20100401 6594 integer :: elfd,iatd,iatom,qvecd 6595 real(dp) :: fac 6596 logical :: iwrite 6597 character(len=500) :: msg 6598 !arrays 6599 real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert) 6600 6601 ! ********************************************************************* 6602 6603 d3cart(1,:,:,:,:,:,:) = reshape(blkval(1,:),shape = (/3,mpert,3,mpert,3,mpert/)) 6604 d3cart(2,:,:,:,:,:,:) = reshape(blkval(2,:),shape = (/3,mpert,3,mpert,3,mpert/)) 6605 6606 !Define a factor to apply if DDB file has been created with the old version of 6607 !the longwave driver. 6608 if (ddb_version <= cvrsio8) then 6609 fac=-two 6610 else 6611 fac=one 6612 end if 6613 6614 !Extraction of quadrupoles (need symmetrization wrt qvecd and elfd) 6615 do iatom = 1,natom 6616 do iatd = 1,3 6617 do elfd = 1,3 6618 do qvecd = 1,elfd-1 6619 if (lwsym==1) then 6620 lwtens(elfd,qvecd,iatd,iatom) = fac * & 6621 (d3cart(2,elfd,natom+2,iatd,iatom,qvecd,natom+8)+d3cart(2,qvecd,natom+2,iatd,iatom,elfd,natom+8)) 6622 lwtens(qvecd,elfd,iatd,iatom) = lwtens(elfd,qvecd,iatd,iatom) 6623 else if (lwsym==0) then 6624 lwtens(elfd,qvecd,iatd,iatom) = fac * d3cart(2,elfd,natom+2,iatd,iatom,qvecd,natom+8) 6625 lwtens(qvecd,elfd,iatd,iatom) = fac * d3cart(2,qvecd,natom+2,iatd,iatom,elfd,natom+8) 6626 end if 6627 end do 6628 if (lwsym==1) then 6629 lwtens(elfd,elfd,iatd,iatom) = fac * two*d3cart(2,elfd,natom+2,iatd,iatom,elfd,natom+8) 6630 else if (lwsym==0) then 6631 lwtens(elfd,elfd,iatd,iatom) = fac * d3cart(2,elfd,natom+2,iatd,iatom,elfd,natom+8) 6632 end if 6633 end do 6634 end do 6635 end do 6636 6637 iwrite = ab_out > 0 6638 6639 if (iwrite) then 6640 if (lwsym==1) then 6641 write(msg,*)' atom dir Qxx Qyy Qzz Qyz Qxz Qxy' 6642 call wrtout([ab_out,std_out],msg) 6643 do iatom= 1, natom 6644 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iatom, 'x',lwtens(1,1,1,iatom),lwtens(2,2,1,iatom),lwtens(3,3,1,iatom), & 6645 & lwtens(2,3,1,iatom),lwtens(1,3,1,iatom),lwtens(1,2,1,iatom) 6646 call wrtout([ab_out,std_out],msg) 6647 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iatom, 'y',lwtens(1,1,2,iatom),lwtens(2,2,2,iatom),lwtens(3,3,2,iatom), & 6648 & lwtens(2,3,2,iatom),lwtens(1,3,2,iatom),lwtens(1,2,2,iatom) 6649 call wrtout([ab_out,std_out],msg) 6650 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iatom, 'z',lwtens(1,1,3,iatom),lwtens(2,2,3,iatom),lwtens(3,3,3,iatom), & 6651 & lwtens(2,3,3,iatom),lwtens(1,3,3,iatom),lwtens(1,2,3,iatom) 6652 call wrtout([ab_out,std_out],msg) 6653 end do 6654 else if (lwsym==0) then 6655 write(msg,*) & 6656 & ' atom dir Pxx Pyy Pzz Pyz Pxz Pxy Pzy Pzx Pyx' 6657 call wrtout([ab_out,std_out],msg) 6658 do iatom= 1, natom 6659 write(msg,'(2x,i3,3x,a3,2x,9f12.6)') iatom, 'x',lwtens(1,1,1,iatom),lwtens(2,2,1,iatom),lwtens(3,3,1,iatom), & 6660 & lwtens(2,3,1,iatom),lwtens(1,3,1,iatom),lwtens(1,2,1,iatom), & 6661 & lwtens(3,2,1,iatom),lwtens(3,1,1,iatom),lwtens(2,1,1,iatom) 6662 call wrtout([ab_out,std_out],msg) 6663 write(msg,'(2x,i3,3x,a3,2x,9f12.6)') iatom, 'y',lwtens(1,1,2,iatom),lwtens(2,2,2,iatom),lwtens(3,3,2,iatom), & 6664 & lwtens(2,3,2,iatom),lwtens(1,3,2,iatom),lwtens(1,2,2,iatom), & 6665 & lwtens(3,2,2,iatom),lwtens(3,1,2,iatom),lwtens(2,1,2,iatom) 6666 call wrtout([ab_out,std_out],msg) 6667 write(msg,'(2x,i3,3x,a3,2x,9f12.6)') iatom, 'z',lwtens(1,1,3,iatom),lwtens(2,2,3,iatom),lwtens(3,3,3,iatom), & 6668 & lwtens(2,3,3,iatom),lwtens(1,3,3,iatom),lwtens(1,2,3,iatom), & 6669 & lwtens(3,2,3,iatom),lwtens(3,1,3,iatom),lwtens(2,1,3,iatom) 6670 call wrtout([ab_out,std_out],msg) 6671 end do 6672 endif 6673 end if 6674 6675 end subroutine dtqdrp
m_ddb/gamma9 [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
gamma9
FUNCTION
This small routine checks if the wavevector qphon and the corresponding normalisation factor represent a phonon at Gamma.
INPUTS
qphon(3)=wavevector qphnrm=normalisation factor qtol=tolerance
OUTPUT
gamma= if 1, means that the wavevector is indeed at Gamma otherwise 0.
SOURCE
1595 subroutine gamma9(gamma,qphon,qphnrm,qtol) 1596 1597 !Arguments ------------------------------- 1598 !scalars 1599 integer,intent(out) :: gamma 1600 real(dp),intent(in) :: qphnrm,qtol 1601 !arrays 1602 real(dp),intent(in) :: qphon(3) 1603 1604 ! ********************************************************************* 1605 1606 if( (abs(qphon(1))<qtol .and. abs(qphon(2))<qtol .and. abs(qphon(3))<qtol) .or. abs(qphnrm)<qtol ) then 1607 gamma=1 1608 else 1609 gamma=0 1610 end if 1611 1612 end subroutine gamma9
m_ddb/lwcart [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
lwcart
FUNCTION
Transform the 3rd-order energy derivative read from the ddb file generated by a long wave calculation into cartesian coordinates, and also...
INPUTS
blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte has been calculated ; 0 otherwise ) d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1) mpert =maximum number of ipert natom= number of atoms rprimd(3,3)=dimensional primitive translations (bohr)
OUTPUT
carflg(3,mpert,3,mpert,3,mpert)=1 if the element of d3cart has been calculated, 0 otherwise d3cart(2,3,mpert,3,mpert,3,mpert)=matrix of third-order energy derivatives in cartesian coordinates
SOURCE
6484 subroutine lwcart(blkflg,carflg,d3,d3cart,gprimd,mpert,natom,rprimd) 6485 6486 !Arguments ------------------------------- 6487 !scalars 6488 integer,intent(in) :: mpert,natom 6489 !arrays 6490 integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert) 6491 integer,intent(out) :: carflg(3,mpert,3,mpert,3,mpert) 6492 real(dp),intent(in) :: d3(2,3,mpert,3,mpert,3,mpert),gprimd(3,3),rprimd(3,3) 6493 real(dp),intent(out) :: d3cart(2,3,mpert,3,mpert,3,mpert) 6494 6495 !Local variables ------------------------- 6496 !scalars 6497 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 6498 integer :: ii 6499 !arrays 6500 integer :: flg1(3),flg2(3) 6501 real(dp) :: vec1(3),vec2(3) 6502 6503 ! ******************************************************************* 6504 6505 !Transform to cartesian coordinates 6506 d3cart(:,:,:,:,:,:,:) = d3(:,:,:,:,:,:,:) 6507 carflg(:,:,:,:,:,:) = 0 6508 6509 do i1pert = 1, mpert 6510 do i2pert = 1, mpert 6511 do i3pert = 1, mpert 6512 6513 do i2dir = 1, 3 6514 do i3dir = 1, 3 6515 do ii= 1, 2 6516 vec1(:) = d3cart(ii,:,i1pert,i2dir,i2pert,i3dir,i3pert) 6517 flg1(:) = blkflg(:,i1pert,i2dir,i2pert,i3dir,i3pert) 6518 call cart39(flg1,flg2,gprimd,i1pert,natom,rprimd,vec1,vec2) 6519 d3cart(ii,:,i1pert,i2dir,i2pert,i3dir,i3pert) = vec2(:) 6520 carflg(:,i1pert,i2dir,i2pert,i3dir,i3pert) = flg2(:) 6521 end do 6522 end do 6523 end do 6524 6525 do i1dir = 1, 3 6526 do i3dir = 1, 3 6527 do ii= 1, 2 6528 vec1(:) = d3cart(ii,i1dir,i1pert,:,i2pert,i3dir,i3pert) 6529 flg1(:) = blkflg(i1dir,i1pert,:,i2pert,i3dir,i3pert) 6530 call cart39(flg1,flg2,gprimd,i2pert,natom,rprimd,vec1,vec2) 6531 d3cart(ii,i1dir,i1pert,:,i2pert,i3dir,i3pert) = vec2(:) 6532 carflg(i1dir,i1pert,:,i2pert,i3dir,i3pert) = flg2(:) 6533 end do 6534 end do 6535 end do 6536 6537 do i1dir = 1, 3 6538 do i2dir = 1, 3 6539 do ii= 1, 2 6540 vec1(:) = d3cart(ii,i1dir,i1pert,i2dir,i2pert,:,i3pert) 6541 flg1(:) = blkflg(i1dir,i1pert,i2dir,i2pert,:,i3pert) 6542 call cart39(flg1,flg2,gprimd,i3pert,natom,rprimd,vec1,vec2) 6543 d3cart(ii,i1dir,i1pert,i2dir,i2pert,:,i3pert) = vec2(:) 6544 carflg(i1dir,i1pert,i2dir,i2pert,:,i3pert) = flg2(:) 6545 end do 6546 end do 6547 end do 6548 6549 end do 6550 end do 6551 end do 6552 6553 end subroutine lwcart
m_ddb/merge_ddb [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
merge_ddb
FUNCTION
Read a list of ddb files and merge them into a single ddb object.
INPUTS
nddb=number of DDBs to merge
filenames=names of input DDB files
outfile=name of the merged DDB file to be written
dscrpt=string description of the final ddb.
chkopt=option for consistency checks between DDB files
(0 --> do not check header consistency between files)
(1 --> check header consistency between files)
OUTPUT
SOURCE
6159 subroutine merge_ddb(nddb, filenames, outfile, dscrpt, chkopt) 6160 6161 !Arguments ------------------------------- 6162 !scalars 6163 integer,intent(in) :: nddb 6164 integer,intent(in) :: chkopt 6165 !arrays 6166 character(len=fnlen),intent(in) :: filenames(nddb) 6167 character(len=fnlen),intent(in) :: outfile, dscrpt 6168 6169 !Local variables ------------------------- 6170 !scalars 6171 integer,parameter :: master=0 6172 integer :: iddb, iddb_mkpt, iddb_psps 6173 integer :: dimekb, matom, mband, mblok, mkpt, nsppol 6174 integer :: mtypat, lmnmax, usepaw, mblktyp, msym 6175 integer :: msize, msize_, mpert 6176 integer :: nblok, iblok, iblok1, iblok2 6177 integer :: comm 6178 logical :: eig2d, can_merge 6179 integer,parameter :: prtvol=0 6180 character(len=500) :: msg 6181 type(ddb_type) :: ddb, ddb2 6182 type(ddb_hdr_type) :: ddb_hdr, ddb_hdr2 6183 type(crystal_t) :: crystal 6184 6185 ! ************************************************************************ 6186 6187 comm = xmpi_world 6188 6189 ! ----------------------------------------------- 6190 ! Read all headers and evaluate arrays dimensions 6191 ! ----------------------------------------------- 6192 if (xmpi_comm_rank(comm) == master) then 6193 call wrtout(std_out, sjoin(ch10, " merge_ddb: Reading all headers.")) 6194 end if 6195 6196 dimekb=0 ; matom=0 ; mband=0 ; mblok=0 ; mkpt=0 ; mpert=0 6197 msize=0 ; mtypat=0 ; lmnmax=0 ; usepaw=0 ; mblktyp=1 6198 iddb_mkpt = 1 ; iddb_psps = nddb 6199 msym=192 6200 6201 eig2d = .False. 6202 do iddb=1,nddb 6203 6204 call ddb_hdr%open_read(filenames(iddb), comm, dimonly=1) 6205 6206 matom=max(matom,ddb_hdr%matom) 6207 6208 ! GA: Should get mkpt from the ddb containing d2eig, if any. 6209 ! In facts, since I removed comparison on the k-points 6210 ! in ddb_hdr_compare, I should add a check to make sure 6211 ! k-points are consistent when merging d2eig data. 6212 if (ddb_hdr%mkpt > mkpt) then 6213 mkpt = ddb_hdr%mkpt 6214 iddb_mkpt = iddb 6215 end if 6216 !mkpt=max(mkpt,ddb_hdr%mkpt) 6217 mtypat=max(mtypat,ddb_hdr%mtypat) 6218 msym=max(msym,ddb_hdr%msym) 6219 mband=max(mband,ddb_hdr%mband) 6220 dimekb=max(dimekb,ddb_hdr%psps%dimekb) 6221 lmnmax=max(lmnmax,ddb_hdr%psps%lmnmax) 6222 usepaw=max(usepaw,ddb_hdr%usepaw) 6223 nsppol = ddb_hdr%nsppol 6224 6225 ! Count the blocks 6226 mblok=mblok+ddb_hdr%nblok 6227 6228 ! Figure out if we are merging eig2d files 6229 if (is_type_d2eig(ddb_hdr%mblktyp)) then 6230 eig2d = .True. 6231 end if 6232 6233 ! Figure out if we are merging d3E blocks and compute msize accordingly 6234 mpert = max(mpert,ddb_hdr%mpert) 6235 msize_ = 3 * mpert * 3 * mpert 6236 if (is_type_d3E(ddb_hdr%mblktyp)) msize_ = msize_ * 3 * mpert 6237 msize = max(msize, msize_) 6238 6239 if (ddb_hdr%with_psps>0 .or. ddb_hdr%psps%usepaw > 0) then 6240 iddb_psps = iddb 6241 end if 6242 6243 end do 6244 6245 ddb%nsppol = nsppol 6246 6247 ! --------------- 6248 ! Allocate arrays 6249 ! --------------- 6250 if (eig2d) then 6251 ! GA: We need to multiply mband by nsppol (to keep array rank below 8) 6252 call ddb%malloc(msize, mblok, matom, mtypat, mpert, mkpt, mband * nsppol) 6253 else 6254 call ddb%malloc(msize, mblok, matom, mtypat, mpert) 6255 end if 6256 6257 ! ------------------------------------------------------- 6258 ! Initialize the output ddb_hdr using the first input ddb 6259 ! ------------------------------------------------------- 6260 6261 ! GA: The last ddb is usually the one that contains the most info on pseudos 6262 ! however, we should check them all and figure out which one has 6263 ! the most info. 6264 6265 call ddb_hdr%free() ! GA: why do I need this? Try to remove 6266 call ddb_hdr%open_read(filenames(1), comm, & 6267 matom=matom,mtypat=mtypat,mband=mband,mkpt=mkpt,& 6268 msym=msym,dimekb=dimekb,lmnmax=lmnmax,usepaw=usepaw) 6269 call ddb_hdr%close() 6270 ddb_hdr%mpert = mpert 6271 ddb_hdr%msize = msize 6272 6273 ! GA: We are setting mkpt at initialization, 6274 ! but netcdf file declares dimension with nkpt. 6275 ! TODO: Should check consistency of nkpt 6276 ! among of all blocks containing eig2d data. 6277 6278 ! ================== 6279 ! Read all databases 6280 ! ================== 6281 6282 nblok = 0 6283 6284 do iddb=1,nddb 6285 6286 ! Open the corresponding input DDB, and read the database file information 6287 write(msg, '(a,a,i6)' )ch10,' read the input derivative database number',iddb 6288 call wrtout(std_out,msg) 6289 6290 ! Note: it is necessary to specify mkpt, otherwise the comparison will crash 6291 call ddb_hdr2%open_read(filenames(iddb), comm, & 6292 matom=matom,mtypat=mtypat,mband=mband,mkpt=mkpt,& 6293 msym=msym,dimekb=dimekb,lmnmax=lmnmax,usepaw=usepaw) 6294 call ddb_hdr2%close() 6295 6296 if (chkopt==1)then 6297 6298 ! Compare the current DDB and input DDB information. 6299 ! In case of an inconsistency, halt the execution. 6300 call wrtout(std_out, ' compare the current and input DDB information') 6301 6302 6303 ! GA: Maybe the problem is that we are comparing uninitialized pawtab 6304 call ddb_hdr%compare(ddb_hdr2) 6305 6306 else 6307 ! No comparison between the current DDB and input DDB information. 6308 call wrtout(std_out,msg) 6309 write(msg, '(3a)' )& 6310 'No comparison/check is performed for the current and input DDB information ',ch10,& 6311 'because argument --nostrict was passed to the command line. ' 6312 ABI_COMMENT(msg) 6313 end if 6314 6315 if (chkopt==1 .or. usepaw==1) then 6316 call ddb_hdr%copy_missing_variables(ddb_hdr2) 6317 end if 6318 6319 ! GA: In principle, this could be done only once, 6320 ! but I could not managed to do that without failing test v8[07]. 6321 if (iddb == iddb_psps) then 6322 call ddb_hdr%copy_psps_from(ddb_hdr2) 6323 end if 6324 6325 call ddb_hdr2%free() 6326 6327 ! Now read the whole DDB 6328 call ddb2%from_file(filenames(iddb), ddb_hdr2, crystal, comm, prtvol, raw=1) 6329 call crystal%free() 6330 call ddb_hdr2%free() 6331 6332 ! -------------------------------------------------------------- 6333 ! Double loop over the blocks of the last ddb and the output ddb 6334 ! -------------------------------------------------------------- 6335 6336 do iblok2=1,ddb2%nblok 6337 6338 can_merge = .false. 6339 do iblok1=1, nblok 6340 6341 can_merge = ddb%can_merge_blocks(ddb2, iblok1, iblok2) 6342 6343 if (can_merge) then 6344 write(msg, '(a,i5,a,a)' )' merge block #',iblok2,' from file ', filenames(iddb) 6345 call wrtout(std_out,msg) 6346 iblok = iblok1 ! Merge with previous block 6347 exit 6348 end if 6349 end do 6350 6351 if (.not. can_merge) then 6352 write(msg, '(a,i5,a,a)' )' add block #',iblok2,' from file ', filenames(iddb) 6353 call wrtout(std_out,msg) 6354 nblok = nblok + 1 6355 iblok = nblok 6356 end if 6357 6358 call ddb%merge_blocks(ddb2, iblok, iblok2) 6359 6360 end do ! iblok2 6361 6362 ! Free memory 6363 call ddb2%free() 6364 6365 end do ! iddb 6366 6367 6368 6369 ddb_hdr%nblok = nblok 6370 ddb%nblok = nblok 6371 ddb_hdr%dscrpt = dscrpt 6372 6373 call ddb_hdr%set_typ(ddb%nblok, ddb%typ) 6374 6375 ddb_hdr%mpert = mpert ! This is done anyway at writing 6376 6377 ! Summarize the merging phase 6378 write(msg, '(a,i6,a)' )' Final DDB has ',nblok,' blocks.' 6379 call wrtout(std_out,msg) 6380 6381 ! Always use format specified with output filename. 6382 ! GA: Might need an extra variable to enforce a different iomode 6383 ddb_hdr%iomode = iomode_from_fname(outfile) 6384 6385 ! GA: This is because netcdf format has more info than txt 6386 ! and psps might be initialized even if with_psps==0 6387 ! Very weird that I have to do this. 6388 ! TODO Do not enforce with_psps=1. Change the test reference instead. 6389 if (ddb_hdr%iomode/=IO_MODE_ETSF) then 6390 if (ddb_hdr%with_psps==0) then 6391 ddb_hdr%psps%dimekb = 0 6392 ddb_hdr%psps%lmnmax = 0 6393 ABI_SFREE(ddb_hdr%psps%ekb) 6394 ABI_SFREE(ddb_hdr%psps%indlmn) 6395 ABI_MALLOC(ddb_hdr%psps%ekb,(ddb_hdr%psps%dimekb,ddb_hdr%mtypat)) 6396 ABI_MALLOC(ddb_hdr%psps%indlmn,(6,ddb_hdr%psps%lmnmax,ddb_hdr%mtypat)) 6397 ddb_hdr%psps%ekb = zero 6398 ddb_hdr%psps%indlmn = zero 6399 end if 6400 end if 6401 6402 ! Enforce full initialization, regardless of DDB content 6403 ddb_hdr%with_psps=1 6404 ddb_hdr%with_dfpt_vars=1 6405 6406 ! Write the final ddb to file. 6407 if (.not. eig2d) then 6408 call ddb%write(ddb_hdr, outfile) 6409 end if 6410 6411 ! ================================= 6412 ! Second derivatives of eigenvalues 6413 ! ================================= 6414 if (eig2d) then 6415 6416 ! GA: Here we assume that the blocks are complete wrt perturbations. 6417 ! No merging of blocks occurs. 6418 ! TODO: Implement merging of partial d2eig blocks 6419 6420 ddb%kpt(:,:) = ddb_hdr%kpt(:,:) 6421 6422 ! Open the output DDB and write the header 6423 call ddb_hdr%open_write(outfile, with_psps=1, comm=comm) 6424 6425 iblok = 0 6426 do iddb=1,nddb 6427 6428 call ddb_hdr2%open_read(filenames(iddb), comm, & 6429 matom=matom,mtypat=mtypat,mband=mband,mkpt=mkpt,& 6430 msym=msym,dimekb=dimekb,lmnmax=lmnmax,usepaw=usepaw) 6431 6432 do iblok2=1,ddb_hdr2%nblok 6433 6434 ! Handle one block at a time 6435 iblok = iblok + 1 6436 call ddb%read_d2eig(ddb_hdr2, iblok, iblok2) 6437 call ddb%write_d2eig(ddb_hdr, iblok) 6438 6439 end do 6440 6441 call ddb_hdr2%close() ! Close the file 6442 call ddb_hdr2%free() ! Free memory 6443 6444 end do 6445 6446 call ddb_hdr%close() 6447 6448 end if 6449 6450 ! ----------- 6451 ! Free memory 6452 ! ----------- 6453 call ddb_hdr%free() 6454 call ddb%free() 6455 6456 end subroutine merge_ddb
m_ddb/nlopt [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
nlopt
FUNCTION
Output of all quantities related to third-order derivatives of the energy. Compute the permutations of the three perturbations, then write out the whole matrix of third order derivatives in reduced coordinates. Finally, compute the non-linear optical susceptibility d and the first-order change in the dielectric susceptibility tensor induced by an atomic displacement.
INPUTS
blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte has been calculated ; 0 otherwise ) d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE gprimd(3,3)=dimensional primitive translations for reciprocal space(bohr^-1) mpert =maximum number of ipert natom= number of atoms rprimd(3,3)=dimensional primitive translations (bohr) ucvol=unit cell volume (bohr^3)
OUTPUT
carflg(3,mpert,3,mpert,3,mpert)=1 if the element of d3cart has been calculated, 0 otherwise d3cart(2,3,mpert,3,mpert,3,mpert)=matrix of third-order energy derivatives in cartesian coordinates
SOURCE
2263 subroutine nlopt(blkflg,carflg,d3,d3cart,gprimd,mpert,natom,rprimd,ucvol) 2264 2265 !Arguments ------------------------------- 2266 !scalars 2267 integer,intent(in) :: mpert,natom 2268 real(dp),intent(in) :: ucvol 2269 !arrays 2270 integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert) 2271 integer,intent(out) :: carflg(3,mpert,3,mpert,3,mpert) 2272 real(dp),intent(in) :: d3(2,3,mpert,3,mpert,3,mpert),gprimd(3,3),rprimd(3,3) 2273 real(dp),intent(out) :: d3cart(2,3,mpert,3,mpert,3,mpert) 2274 2275 !Local variables ------------------------- 2276 !scalars 2277 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 2278 !arrays 2279 integer :: flg1(3),flg2(3) 2280 real(dp) :: vec1(3),vec2(3) 2281 2282 ! ******************************************************************* 2283 2284 !Compute the permutations of the perturbations 2285 2286 d3cart(:,:,:,:,:,:,:) = 0._dp 2287 2288 do i1pert = 1,mpert 2289 do i2pert = 1,mpert 2290 do i3pert = 1,mpert 2291 do i1dir=1,3 2292 do i2dir=1,3 2293 do i3dir=1,3 2294 2295 ! Check if all elements are available 2296 2297 if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=0).and. & 2298 (blkflg(i1dir,i1pert,i3dir,i3pert,i2dir,i2pert)/=0).and. & 2299 (blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert)/=0).and. & 2300 (blkflg(i2dir,i2pert,i3dir,i3pert,i1dir,i1pert)/=0).and. & 2301 (blkflg(i3dir,i3pert,i1dir,i1pert,i2dir,i2pert)/=0).and. & 2302 (blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert)/=0)) then 2303 2304 d3cart(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = & 2305 ( d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) + & 2306 d3(:,i1dir,i1pert,i3dir,i3pert,i2dir,i2pert) + & 2307 d3(:,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) + & 2308 d3(:,i2dir,i2pert,i3dir,i3pert,i1dir,i1pert) + & 2309 d3(:,i3dir,i3pert,i1dir,i1pert,i2dir,i2pert) + & 2310 d3(:,i3dir,i3pert,i2dir,i2pert,i1dir,i1pert))*sixth 2311 2312 end if 2313 end do 2314 end do 2315 end do 2316 end do 2317 end do 2318 end do 2319 2320 !Transform to cartesian coordinates 2321 carflg(:,:,:,:,:,:) = 0 2322 2323 do i1pert = 1, mpert 2324 do i2pert = 1, mpert 2325 do i3pert = 1, mpert 2326 2327 do i2dir = 1, 3 2328 do i3dir = 1, 3 2329 2330 vec1(:) = d3cart(1,:,i1pert,i2dir,i2pert,i3dir,i3pert) 2331 flg1(:) = blkflg(:,i1pert,i2dir,i2pert,i3dir,i3pert) 2332 call cart39(flg1,flg2,gprimd,i1pert,natom,rprimd,vec1,vec2) 2333 d3cart(1,:,i1pert,i2dir,i2pert,i3dir,i3pert) = vec2(:) 2334 carflg(:,i1pert,i2dir,i2pert,i3dir,i3pert) = flg2(:) 2335 2336 end do 2337 end do 2338 2339 do i1dir = 1, 3 2340 do i3dir = 1, 3 2341 vec1(:) = d3cart(1,i1dir,i1pert,:,i2pert,i3dir,i3pert) 2342 flg1(:) = blkflg(i1dir,i1pert,:,i2pert,i3dir,i3pert) 2343 call cart39(flg1,flg2,gprimd,i2pert,natom,rprimd,vec1,vec2) 2344 d3cart(1,i1dir,i1pert,:,i2pert,i3dir,i3pert) = vec2(:) 2345 carflg(i1dir,i1pert,:,i2pert,i3dir,i3pert) = flg2(:) 2346 end do 2347 end do 2348 2349 do i1dir = 1, 3 2350 do i2dir = 1, 3 2351 vec1(:) = d3cart(1,i1dir,i1pert,i2dir,i2pert,:,i3pert) 2352 flg1(:) = blkflg(i1dir,i1pert,i2dir,i2pert,:,i3pert) 2353 call cart39(flg1,flg2,gprimd,i3pert,natom,rprimd,vec1,vec2) 2354 d3cart(1,i1dir,i1pert,i2dir,i2pert,:,i3pert) = vec2(:) 2355 carflg(i1dir,i1pert,i2dir,i2pert,:,i3pert) = flg2(:) 2356 end do 2357 end do 2358 2359 end do 2360 end do 2361 end do 2362 2363 ! Compute non linear-optical coefficients d_ijk (atomic units) 2364 i1pert = natom+2 2365 d3cart(:,:,i1pert,:,i1pert,:,i1pert) = -3._dp*d3cart(:,:,i1pert,:,i1pert,:,i1pert)/(ucvol*2._dp) 2366 2367 ! Compute first-order change in the electronic dielectric 2368 ! susceptibility (Bohr^-1) induced by an atomic displacement 2369 d3cart(1:2,1:3,1:natom,1:3,natom + 2,1:3,natom + 2) = -6._dp*d3cart(1:2,1:3,1:natom,1:3,natom + 2,1:3,natom + 2)/ucvol 2370 2371 end subroutine nlopt
m_ddb/rdddb9 [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
rdddb9
FUNCTION
This routine reads the derivative database entirely, for use in ppddb9, and performs some checks and symmetrisation At the end, the whole DDB is in central memory, contained in the array ddb%val(2,msize,ddb%nblok). The information on it is contained in the four arrays ddb%flg(msize,ddb%nblok) : blok flag for each element ddb%qpt(9,ddb%nblok) : blok wavevector (unnormalized) ddb%nrm(3,ddb%nblok) : blok wavevector normalization ddb%typ(ddb%nblok) : blok type
INPUTS
unddb = unit number for DDB io
dimekb=dimension of ekb (for the time being, only for norm- conserving psps)
iout=unit number for output of formatted data
lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
=if useylm=0, max number of (l,n) comp. over all type of psps
mband=maximum number of bands
mpert =maximum number of ipert
msize=maximum size of data blocks
msym =maximum number of symmetry elements in space group
natom = number of atoms
ntypat=number of atom types
usepaw= 0 for non paw calculation; =1 for paw calculation
[raw] = 1 -> do not perform any symetrization or transformation to cartesian coordinates.
0 (default) -> do perform these transformations.
OUTPUT
acell(3)=length scales of cell (bohr)
amu(ntypat)=mass of the atoms (atomic mass unit)
ddb: ddb blok datatype
contents: ddb%flg(msize,nblok)= flag of existence for each element of the DDB
ddb%nrm(3,nblok) : blok wavevector normalization
ddb%qpt(9,nblok) : blok wavevector (unnormalized)
ddb%typ(nblok) : blok type
ddb%val(2,msize,nblok)= value of each complex element of the DDB
ddb%nblok= number of bloks in the DDB
gmet(3,3)=reciprocal space metric tensor in bohr**-2
gprim(3,3)=dimensionless reciprocal space primitive translations
indsym(4,msym,natom)=indirect indexing array for symmetries
natom=number of atoms in cell
nsym=number of space group symmetries
rmet(3,3)=metric tensor in real space (bohr^2)
rprim(3,3)= primitive translation vectors
symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
symafm(nsym)=Anti-ferromagnetic symmetries.
tnons(3,nsym)=fractional nonsymmorphic translations
typat(natom)=type integer for each atom in cell
ucvol=unit cell volume in bohr**3
xcart(3,natom)=atomic cartesian coordinates
xred(3,natom)=fractional dimensionless atomic coordinates
zion(ntypat)=charge on each type of atom (real number)
znucl(ntypat)=Nuclear charge for each type of pseudopotential
SOURCE
2048 subroutine rdddb9(ddb,ddb_hdr,unddb,& 2049 acell,amu,gmet,gprim,indsym,& 2050 mband,mpert,msize,msym,natom,nkpt,nsym,ntypat,& 2051 rmet,rprim,symrec,symrel,symafm,tnons,typat,ucvol,& 2052 xcart,xred,zion,znucl,raw) 2053 2054 !Arguments ------------------------------- 2055 ! NOTE: these are used for dimensioning and then re-assigned in ioddb8. 2056 ! This is almost definitely bad practice. In particular 2057 ! it should be indsym(4,msym,natom), 2058 ! and 2059 ! the allocation allocate(kpt(3,nkpt)) is strange 2060 !scalars 2061 integer,intent(in) :: unddb,mband,mpert,msize,msym 2062 integer,intent(inout) :: natom,nkpt,nsym,ntypat 2063 real(dp),intent(out) :: ucvol 2064 type(ddb_type),intent(inout) :: ddb 2065 type(ddb_hdr_type),intent(inout) :: ddb_hdr 2066 integer,optional,intent(in) :: raw 2067 !arrays 2068 integer,intent(inout) :: indsym(4,msym,natom) 2069 integer,intent(out) :: symrec(3,3,msym),symrel(3,3,msym),symafm(msym) 2070 integer,intent(out) :: typat(natom) 2071 real(dp),intent(out) :: acell(3),amu(ntypat) 2072 real(dp),intent(out) :: gmet(3,3),gprim(3,3),rmet(3,3) 2073 real(dp),intent(out) :: rprim(3,3),tnons(3,msym),xcart(3,natom),xred(3,natom) 2074 real(dp),intent(out) :: zion(ntypat),znucl(ntypat) 2075 2076 !Local variables ------------------------- 2077 !mtyplo=maximum number of type, locally 2078 !scalars 2079 integer,parameter :: msppol=2,mtyplo=6 2080 integer :: raw_ 2081 integer :: iblok,isym 2082 real(dp),parameter :: tolsym8=tol8 2083 !arrays 2084 real(dp) :: gprimd(3,3),rprimd(3,3) 2085 2086 ! ********************************************************************* 2087 2088 DBG_ENTER("COLL") 2089 2090 if (present(raw)) then 2091 raw_ = raw 2092 else 2093 raw_ = 0 2094 end if 2095 2096 ! FIXME 2097 ! GA: Most of this stuff could be moved up to the calling routine 2098 2099 nsym = ddb_hdr%nsym 2100 acell = ddb_hdr%acell 2101 rprim = ddb_hdr%rprim 2102 2103 amu(:) = ddb_hdr%amu(1:ntypat) 2104 typat(:) = ddb_hdr%typat(1:natom) 2105 zion(:) = ddb_hdr%zion(1:ntypat) 2106 znucl(:) = ddb_hdr%znucl(1:ntypat) 2107 2108 symafm(:) = ddb_hdr%symafm(:) 2109 symrel(:,:,:) = ddb_hdr%symrel(:,:,:) 2110 tnons(:,:) = ddb_hdr%tnons(:,:) 2111 2112 xred(:,:) = ddb_hdr%xred(:,:) 2113 2114 !call ddb_hdr%free() 2115 2116 ! Compute different matrices in real and reciprocal space, also 2117 ! checks whether ucvol is positive. 2118 call mkrdim(acell,rprim,rprimd) 2119 2120 ! call metric without printing to output 2121 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2122 2123 ! Obtain reciprocal space primitive transl g from inverse trans of r 2124 ! (Unlike in abinit, gprim is used throughout ifc; should be changed, later) 2125 call matr3inv(rprim,gprim) 2126 2127 ! Generate atom positions in cartesian coordinates 2128 call xred2xcart(natom,rprimd,xcart,xred) 2129 2130 ! Transposed inversion of the symmetry matrices, for use in the reciprocal space 2131 do isym=1,nsym 2132 call mati3inv(symrel(:,:,isym),symrec(:,:,isym)) 2133 end do 2134 2135 ! SYMATM generates for all the atoms and all the symmetries, the atom 2136 ! on which the referenced one is sent and also the translation bringing 2137 ! back this atom to the referenced unit cell 2138 ! GA: symatm was already called in crystal_init, no need to do it again. 2139 call symatm(indsym,natom,nsym,symrec,tnons,tolsym8,typat,xred) 2140 2141 !write(msg, '(3a,i0,a)' )ch10,ch10,' rdddb9: read ',ddb%nblok,' blocks from the input DDB ' 2142 !call wrtout(std_out,msg) 2143 2144 ! Read the blocks from the input database, and close it. 2145 do iblok=1,ddb%nblok 2146 2147 call ddb%read_block_txt(iblok,mband,mpert,msize,nkpt,unddb) 2148 2149 if (raw_ == 0) then 2150 call ddb%symmetrize_and_transform(ddb_hdr%crystal,iblok) 2151 end if 2152 2153 end do ! iblok 2154 2155 DBG_EXIT("COLL") 2156 2157 end subroutine rdddb9
m_ddb/symdm9 [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
symdm9
FUNCTION
Use the set of special k points calculated by the Monkhorst & Pack Technique. Check if all the information for the k points are present in the DDB to determine their dynamical matrices. Generate the dynamical matrices of the set of k points which samples homogeneously the entire Brillouin zone.
INPUTS
%flg(nsize,nblok)= flag of existence for each element of the DDB %nrm(1,nblok)=norm of qpt providing normalization %qpt(1<ii<9,nblok)=q vector of a phonon mode (ii=1,2,3) %typ(nblok)=1 or 2 depending on non-stationary or stationary block 3 for third order derivatives %val(2,3*mpert*3*mpert,nblok)= all the dynamical matrices gprimd(3,3)=dimensionlal primitive translations in reciprocal space indsym = mapping of atoms under symops mpert =maximum number of ipert natom=number of atoms in unit cell %nblok=number of blocks in the DDB nqpt=number of special q points nsym=number of space group symmetries rfmeth = 1 or -1 if non-stationary block 2 or -2 if stationary block 3 or -3 if third order derivatives 85 if molecular Berry curvature positive if symmetries are used to set elements to zero whenever possible, negative to prevent this to happen. rprimd(3,3)=dimensional primitive translations in real space spqpt(3,nqpt)=set of special q points generated by the Monkhorst & Pack Method symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space) comm=MPI communicator.
OUTPUT
dynmat(2,3,natom,3,natom,nqpt)=dynamical matrices relative to the q points of the B.Z. sampling
[qmissing]=Allocatable array with the indices of the q-points in the BZ that could not be obtained
by symmetry. If qmissing is present, the routine does not stop if the full BZ cannot be reconstructed.
The caller is responsible for filling the missing entries.
TODO
* A full description of the inputs should be included
NOTES
Time-reversal symmetry is always assumed The time-reversal is correctly used for the MBC (Molecular Berry curvature): G(-k)=G^*(k)
SOURCE
6804 subroutine symdm9(ddb, dynmat, gprimd, indsym, mpert, natom, nqpt, nsym, rfmeth,& 6805 rprimd, spqpt, symrec, symrel, comm, qmissing) 6806 6807 !Arguments ------------------------------- 6808 !scalars 6809 type(ddb_type),intent(in) :: ddb 6810 integer,intent(in) :: mpert,natom,nqpt,nsym,rfmeth,comm 6811 !arrays 6812 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 6813 integer,allocatable,optional,intent(out) :: qmissing(:) 6814 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 6815 real(dp),intent(in) :: spqpt(3,nqpt) 6816 real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt) 6817 6818 !Local variables ------------------------- 6819 !scalars 6820 integer :: ia,ib,iblok,idir1,idir2,ii,ipert1,ipert2,iqpt,isym,jj,kk,ll 6821 integer :: mu,nu,q1,q2,nqmiss,nprocs,my_rank,ierr,index 6822 real(dp),parameter :: tol=2.d-8 6823 real(dp) :: sign1, sign2 6824 !tolerance for equality of q points between those of the DDB and those of the sampling grid 6825 real(dp) :: arg1,arg2,im,re,sumi,sumr 6826 logical :: allow_qmiss 6827 character(len=500) :: msg 6828 !arrays 6829 integer,allocatable :: qtest(:,:) 6830 integer :: qmiss_(nqpt) 6831 real(dp) :: qq(3),qsym(6),ss(3,3) 6832 real(dp),allocatable :: ddd(:,:,:,:,:) 6833 6834 ! ********************************************************************* 6835 6836 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 6837 6838 ! Initialize output (some q-points might not be reconstructed if qmissing is present) 6839 dynmat = zero 6840 allow_qmiss = (present(qmissing)) 6841 6842 ABI_MALLOC(ddd,(2,3,natom,3,natom)) 6843 ! Check if the blkqpt points and their symmetrics are sufficient 6844 ! in the DDB to retrieve all the q points of the B.Z. sampling 6845 6846 !Initialization of a test variable 6847 ! qtest(iqpt,1)=iblok 6848 ! qtest(iqpt,2)=isym 6849 ! qtest(iqpt,3)=time_reversal 6850 ABI_MALLOC(qtest,(nqpt,3)) 6851 do iqpt=1,nqpt 6852 qtest(iqpt,1)=0 6853 end do 6854 6855 !Q points coming from the DDB 6856 !write(std_out,*)' Nbr. of Blocks -> ',nblok 6857 ! TODO: This part scales badly with nblock/nqpt 6858 ! One could use listkk or rearrange the loop so that iqpt comes first and then MPI-parallelize. 6859 6860 do iblok=1,ddb%nblok 6861 6862 if (abs(ddb%typ(iblok)) == abs(rfmeth)) then 6863 qq(1)=ddb%qpt(1,iblok)/ddb%nrm(1,iblok) 6864 qq(2)=ddb%qpt(2,iblok)/ddb%nrm(1,iblok) 6865 qq(3)=ddb%qpt(3,iblok)/ddb%nrm(1,iblok) 6866 6867 ! Calculation of the symmetric points (including Time Reversal) 6868 do isym=1,nsym 6869 qsym(1)=qq(1)*symrec(1,1,isym)+qq(2)*symrec(1,2,isym)+qq(3)*symrec(1,3,isym) 6870 qsym(2)=qq(1)*symrec(2,1,isym)+qq(2)*symrec(2,2,isym)+qq(3)*symrec(2,3,isym) 6871 qsym(3)=qq(1)*symrec(3,1,isym)+qq(2)*symrec(3,2,isym)+qq(3)*symrec(3,3,isym) 6872 6873 ! Dont forget the Time Reversal symmetry 6874 qsym(4)=-qq(1)*symrec(1,1,isym)-qq(2)*symrec(1,2,isym)-qq(3)*symrec(1,3,isym) 6875 qsym(5)=-qq(1)*symrec(2,1,isym)-qq(2)*symrec(2,2,isym)-qq(3)*symrec(2,3,isym) 6876 qsym(6)=-qq(1)*symrec(3,1,isym)-qq(2)*symrec(3,2,isym)-qq(3)*symrec(3,3,isym) 6877 6878 ! Comparison between the q points and their symmetric points 6879 ! and the set of q points which samples the entire Brillouin zone 6880 do iqpt=1,nqpt 6881 6882 if (mod(abs(spqpt(1,iqpt)-qsym(1))+tol,1._dp)<2*tol)then 6883 if (mod(abs(spqpt(2,iqpt)-qsym(2))+tol,1._dp)<2*tol)then 6884 if (mod(abs(spqpt(3,iqpt)-qsym(3))+tol,1._dp)<2*tol)then 6885 6886 ! write(std_out,*)' q point from the DDB ! ' 6887 ! write(std_out,*)' block -> ',iblok 6888 ! write(std_out,*)' sym. -> ',isym 6889 ! write(std_out,*)' No Time Reversal ' 6890 ! write(std_out,*)'(',qsym(1),',',qsym(2),',',qsym(3),')' 6891 ! write(std_out,*)' ' 6892 qtest(iqpt,1)=iblok 6893 qtest(iqpt,2)=isym 6894 qtest(iqpt,3)=0 6895 end if 6896 end if 6897 end if 6898 6899 if (mod(abs(spqpt(1,iqpt)-qsym(4))+tol,1._dp)<2*tol)then 6900 if (mod(abs(spqpt(2,iqpt)-qsym(5))+tol,1._dp)<2*tol)then 6901 if (mod(abs(spqpt(3,iqpt)-qsym(6))+tol,1._dp)<2*tol)then 6902 6903 ! write(std_out,*)' q point from the DDB ! ' 6904 ! write(std_out,*)' block -> ',iblok 6905 ! write(std_out,*)' sym. -> ',isym 6906 ! write(std_out,*)' Time Reversal ' 6907 ! write(std_out,*)'(',qsym(4),',',qsym(5),',',qsym(6),')' 6908 ! write(std_out,*)' ' 6909 6910 qtest(iqpt,1)=iblok 6911 qtest(iqpt,2)=isym 6912 qtest(iqpt,3)=1 6913 end if 6914 end if 6915 end if 6916 6917 end do ! iqpt 6918 end do ! isym 6919 6920 end if 6921 end do ! iblok 6922 6923 ! Check if all the information relatives to the q points sampling are found in the DDB; 6924 ! if not => stop message 6925 nqmiss = 0 6926 do iqpt=1,nqpt 6927 if (qtest(iqpt,1)==0) then 6928 nqmiss = nqmiss + 1 6929 qmiss_(nqmiss) = iqpt 6930 write(msg, '(3a)' )& 6931 ' symdm9: the bloks found in the DDB are characterized',ch10,& 6932 ' by the following wavevectors :' 6933 call wrtout(std_out,msg) 6934 do iblok=1,ddb%nblok 6935 write(msg, '(a,4d20.12)')' ',ddb%qpt(1,iblok),ddb%qpt(2,iblok),ddb%qpt(3,iblok),ddb%nrm(1,iblok) 6936 call wrtout(std_out,msg) 6937 end do 6938 write(msg, '(a,a,a,i0,a,a,a,3es16.6,a,a,a,a)' )& 6939 'Information is missing in the DDB.',ch10,& 6940 'The dynamical matrix number ',iqpt,' cannot be built,',ch10,& 6941 'since no block with qpt:',spqpt(1:3,iqpt),ch10,& 6942 'has been found.',ch10,& 6943 'Action: add the required block in the DDB, or modify the q-mesh your input file.' 6944 if (.not.allow_qmiss) then 6945 ABI_ERROR(msg) 6946 else 6947 !continue 6948 ABI_COMMENT(msg) 6949 end if 6950 end if 6951 end do 6952 6953 ! Will return a list with the index of the q-points that could not be symmetrized. 6954 if (allow_qmiss) then 6955 ABI_MALLOC(qmissing, (nqmiss)) 6956 if (nqmiss > 0) qmissing = qmiss_(1:nqmiss) 6957 end if 6958 6959 ! Generation of the dynamical matrices relative to the q points 6960 ! of the set which samples the entire Brillouin zone 6961 do iqpt=1,nqpt 6962 if (mod(iqpt, nprocs) /= my_rank) cycle ! mpi-parallelism 6963 6964 q1=qtest(iqpt,1) ! iblok 6965 q2=qtest(iqpt,2) ! isym 6966 ! Skip this q-point if don't have enough info and allow_qmiss 6967 if (allow_qmiss .and. q1==0) cycle 6968 6969 ! Check if the symmetry accompagnied with time reversal : q <- -q 6970 do ii=1,3 6971 qq(ii)=ddb%qpt(ii,q1)/ddb%nrm(1,q1) 6972 end do 6973 if (qtest(iqpt,3)/=0) qq(:) = -qq(:) 6974 ! 6975 do ii=1,3 6976 do jj=1,3 6977 ss(ii,jj)=zero 6978 do kk=1,3 6979 do ll=1,3 6980 ss(ii,jj) = ss(ii,jj) + rprimd(ii,kk) * symrel(kk,ll,q2) * gprimd(jj,ll) 6981 end do 6982 end do 6983 end do 6984 end do 6985 6986 ! Check whether all the information is contained in the DDB 6987 do ipert2=1,natom 6988 do idir2=1,3 6989 do ipert1=1,natom 6990 do idir1=1,3 6991 index = idir1+ 3*((ipert1-1)+ddb%mpert*((idir2-1)+3*(ipert2-1))) 6992 !if(ddb%flg(idir1,ipert1,idir2,ipert2,q1)/=1)then 6993 if(ddb%flg(index,q1)/=1)then 6994 write(msg, '(a,a,a,i0,a,a,a,4(i0,1x),a,a,a,a)' )& 6995 'Elements are missing in the DDB.',ch10,& 6996 'In block iq1: ',q1,' the following element is missing: ',ch10,& 6997 '(idir1, ipert1, idir2, ipert2): ',idir1,ipert1,idir2,ipert2,ch10,& 6998 'Action: add the required information in the DDB with mrgddb,',ch10,& 6999 'and/or check that all irreducible perturbations have been computed.' 7000 ABI_ERROR(msg) 7001 end if 7002 end do 7003 end do 7004 end do 7005 end do 7006 7007 ! Read the dynamical matrices in the DDB 7008 do ipert2=1,natom 7009 do idir2=1,3 7010 do ipert1=1,natom 7011 do idir1=1,3 7012 ddd(:,idir1,ipert1,idir2,ipert2)=ddb%val(:,idir1+3*(ipert1-1+mpert*(idir2-1+3*(ipert2-1))),q1) 7013 end do 7014 end do 7015 end do 7016 end do 7017 7018 ! determine sign of complex conjugation 7019 ! If there is Time Reversal : D.M. <- Complex Conjugate D.M. 7020 ! MBC <- Complex Conjugate MBC 7021 if (qtest(iqpt,3)==0) then 7022 sign1 = one 7023 sign2 = one 7024 else ! if timrev -> complex conjugate 7025 sign1 = one 7026 sign2 = -one 7027 end if 7028 7029 ! Calculation of the dynamical matrix of a symmetrical q point 7030 do ia=1,natom 7031 do ib=1,natom 7032 ! write(std_out,*)'atom-> ',ia,indsym(4,q2,ia); write(std_out,*)'atom-> ',ib,indsym(4,q2,ib) 7033 arg1=two_pi*(qq(1)*indsym(1,q2,ia)+qq(2)*indsym(2,q2,ia)+qq(3)*indsym(3,q2,ia)) 7034 arg2=two_pi*(qq(1)*indsym(1,q2,ib)+qq(2)*indsym(2,q2,ib)+qq(3)*indsym(3,q2,ib)) 7035 re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2) 7036 im=cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2) 7037 do mu=1,3 7038 do nu=1,3 7039 sumr=zero 7040 sumi=zero 7041 do ii=1,3 7042 do jj=1,3 7043 sumr=sumr+ss(mu,ii)*ss(nu,jj)*ddd(1,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib)) 7044 sumi=sumi+ss(mu,ii)*ss(nu,jj)*ddd(2,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib)) 7045 end do 7046 end do 7047 7048 ! Dynmat -> Dynamical Matrix for the q point of the sampling 7049 ! write(std_out,*)' Sumr -> ',mu,nu,sumr; write(std_out,*)' Sumi -> ',mu,nu,sumi 7050 dynmat(1,mu,ia,nu,ib,iqpt) = sign1*re*sumr - sign2*im*sumi 7051 dynmat(2,mu,ia,nu,ib,iqpt) = sign2*re*sumi + sign1*im*sumr 7052 end do ! coordinates 7053 end do 7054 7055 end do ! ia atoms 7056 end do ! ib atoms 7057 end do ! q points of the sampling 7058 7059 ABI_FREE(ddd) 7060 ABI_FREE(qtest) 7061 7062 7063 call xmpi_sum(dynmat, comm, ierr) 7064 7065 end subroutine symdm9