TABLE OF CONTENTS
- ABINIT/m_ddb_flexo
- ABINIT/m_ddb_flexo/ddb_flexo
- m_ddb/dm_psinv
- m_ddb/dtciflexo
- m_ddb/dtlattflexo
- m_ddb/dtmixflexo
ABINIT/m_ddb_flexo [ Modules ]
NAME
m_ddb_flexo
FUNCTION
FIXME: add description.
COPYRIGHT
Copyright (C) 2019-2024 ABINIT group (MR,MS) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
SOURCE
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 module m_ddb_flexo 25 26 use defs_basis 27 use m_abicore 28 use m_profiling_abi 29 use m_errors 30 31 use m_fstrings, only : itoa,sjoin 32 use m_ddb_hdr 33 use m_ddb 34 use m_crystal, only : crystal_t 35 use m_dynmat, only : asria_corr,cart39 36 37 implicit none 38 39 private 40 41 public :: ddb_flexo 42 ! ************************************************************************* 43 44 contains
ABINIT/m_ddb_flexo/ddb_flexo [ Functions ]
NAME
ddb_flexo
FUNCTION
Get all the contributions to the flexoelectric tensor
INPUTS
asr= if /=0 acustic sume rule is imposed on the dynamical matrix d2asr(2,3,natom,3,natom)=ASR-correction ddb<type(ddb_type)>=2nd order derivative database. ddb<type(ddb_type)>=Long wave 3rd order derivative database. ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files. Crystal<type(crystal_t)>=Crystal structure parameters filnamddb = name of the ddb file flexoflg= 1 -> Computes all contributions to FxE 2 -> Computes electronic (clamped ion) contribution to FxE 3 -> Computes mixed contribution to FxE 4 -> Computes lattice contribution to FxE prtvol= if > 1 print all individual quantities of the lattice contribution FxE zeff(3,3,natom)= Born Effective charges
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
77 subroutine ddb_flexo(asr,d2asr,ddb,ddb_lw,ddb_version,crystal,filnamddb,flexoflg,prtvol,zeff) 78 79 implicit none 80 81 !Arguments ------------------------------------ 82 !scalars 83 integer,intent(in) :: ddb_version 84 integer , intent(in) :: asr,flexoflg,prtvol 85 class(ddb_type),intent(in) :: ddb,ddb_lw 86 type(crystal_t),intent(in) :: crystal 87 character(len=fnlen) :: filnamddb 88 !arrays 89 real(dp),intent(in) :: d2asr(2,3,ddb%natom,3,ddb%natom) 90 real(dp),intent(in) :: zeff(3,3,ddb%natom) 91 92 !Local variables------------------------------- 93 integer :: elfd,iblok,ivar,jblok,kblok,lblok,lwsym,qvecd 94 logical :: intstrn_only,iwrite 95 character(len=500) :: msg 96 97 !arrays 98 integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/) 99 integer :: rfelfd(4),rfphon(4),rfstrs(4) 100 integer :: rfqvec(4) 101 real(dp) :: qphnrm(3),qphon(3,3) 102 real(dp) :: ciflexo(3,3,3,3) 103 real(dp) :: intstrn(3,3,3,ddb%natom) 104 real(dp) :: piezofr(3,ddb%natom,3,3) 105 real(dp) :: lattflexo(3,3,3,3) 106 real(dp) :: mixflexo(3,3,3,3) 107 real(dp) :: pol1(3,3,3,ddb%natom) 108 real(dp) :: psinvdm(3*ddb%natom,3*ddb%natom) 109 real(dp) :: totflexo(3,3,3,3) 110 character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/) 111 112 ! ************************************************************************* 113 114 DBG_ENTER("COLL") 115 116 ! First get the clamped-ion flexoelectric tensor 117 ciflexo(:,:,:,:)=zero 118 if (flexoflg==1.or.flexoflg==2) then 119 120 rfphon(:)=0 121 rfelfd(:)=0 122 rfstrs(:)=0 123 rfqvec(:)=0 124 125 ! Look for the Gamma Block in the DDB 126 qphon(:,:)=zero 127 qphnrm(:)=one 128 rfphon(:)=0 129 rfelfd(1)=2 130 rfstrs(2)=3 131 rfqvec(3)=1 132 133 write(msg, '(2a)' ) ch10," Extract the electronic flexoelectric coeficients from 3DTE" 134 call wrtout(std_out,msg,'COLL') 135 call ddb_lw%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec) 136 137 if (iblok == 0) then 138 call wrtout(std_out, " ") 139 call wrtout(std_out, "--- !WARNING") 140 call wrtout(std_out, sjoin("- Cannot find clamped ion FxE tensor in DDB file:", filnamddb)) 141 call wrtout(std_out, " flexoflag=1 or 2 requires the DDB file to include the corresponding long wave 3rd derivatives") 142 else 143 call dtciflexo(ddb_lw%val(:,:,iblok),ddb_version,ddb%mpert,ddb%natom,ciflexo,crystal%ucvol) 144 end if 145 146 end if 147 148 ! Then get the mixed contribution to the flexoelectric tensor 149 !Activate the calculation of internal strain necessary for lattice mediated contribution. 150 intstrn_only=.false.;if (flexoflg==4) intstrn_only=.true. 151 mixflexo(:,:,:,:)=zero 152 if (flexoflg==1.or.flexoflg==3.or.intstrn_only) then 153 154 ! Extract the P^(1) tensor from the DDB 155 if (.not.intstrn_only) then 156 lwsym=0 157 iblok = ddb_lw%get_quadrupoles(ddb_version,lwsym,BLKTYP_d3E_lw,pol1) 158 end if 159 160 rfphon(:)=0 161 rfelfd(:)=0 162 rfstrs(:)=0 163 rfqvec(:)=0 164 165 ! Look for the Gamma Block of the Phi^(1) tensor in the DDB 166 qphon(:,:)=zero 167 qphnrm(:)=one 168 rfphon(1)=1 169 rfphon(2)=1 170 rfqvec(3)=1 171 172 write(msg, '(2a)' ) ch10," Extract the Phi^(1) coeficients from 3DTE" 173 call wrtout(std_out,msg,'COLL') 174 call ddb_lw%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec) 175 176 if (iblok == 0) then 177 call wrtout(std_out, " ") 178 call wrtout(std_out, "--- !WARNING") 179 call wrtout(std_out, sjoin("- Cannot find Phi^(1) tensor in DDB file:", filnamddb)) 180 call wrtout(std_out, " flexoflag=1 or 3 requires the DDB file to include the corresponding long wave 3rd derivatives") 181 end if 182 183 ! Look for th block that contains the forces 184 qphon(:,:)=zero 185 qphnrm(:)=one 186 rfphon(:)=0 187 rfphon(4)=1 188 189 write(msg, '(2a)' ) ch10," Extract the forces from 1DTE" 190 call wrtout(std_out,msg,'COLL') 191 call ddb%get_block(kblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,4) 192 193 if (kblok == 0) then 194 call wrtout(std_out, " ") 195 call wrtout(std_out, "--- !WARNING") 196 call wrtout(std_out, sjoin("- Cannot find forces in DDB file:", filnamddb)) 197 call wrtout(std_out, " If there are nonzero residual atomic forces on the structure") 198 call wrtout(std_out, " flexoflag=1 or 3 will produce an improper piezoelectric force response tensor") 199 call wrtout(std_out, " and a wrong value for the mixed contribution to the flexoelectric tensor") 200 end if 201 202 ! Look for the Gamma Block of the dynamical matrix in the DDB 203 qphon(:,:)=zero 204 qphnrm(:)=one 205 rfphon(:)=0 206 rfphon(1:2)=1 207 208 write(msg, '(2a)' ) ch10," Extract the Dynamical Matrix from 2DTE" 209 call wrtout(std_out,msg,'COLL') 210 call ddb%get_block(jblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,1) 211 212 if (jblok == 0) then 213 call wrtout(std_out, " ") 214 call wrtout(std_out, "--- !WARNING") 215 call wrtout(std_out, sjoin("- Cannot find Gamma point Dynamical Matrix in DDB file:", filnamddb)) 216 call wrtout(std_out, " flexoflag=1 or 3 requires the DDB file to include the corresponding 2nd derivatives") 217 end if 218 219 if (iblok/=0.and.jblok/=0) then 220 call dtmixflexo(asr,d2asr,ddb%val(:,:,kblok),ddb%val(:,:,jblok),ddb_lw%val(:,:,iblok),ddb_version,crystal%gprimd,& 221 & intstrn,intstrn_only,mixflexo,ddb%mpert,ddb%natom,piezofr,pol1,psinvdm,crystal%rprimd,crystal%ucvol) 222 end if 223 224 end if 225 226 ! Finally get the lattice mediated contribution to the flexoelectric tensor 227 lattflexo(:,:,:,:)=zero 228 if (flexoflg==1.or.flexoflg==4) then 229 230 rfphon(:)=0 231 rfelfd(:)=0 232 rfstrs(:)=0 233 rfqvec(:)=0 234 235 ! Look for the Gamma Block of the Phi^(1) tensor in the DDB 236 qphon(:,:)=zero 237 qphnrm(:)=one 238 rfphon(1)=1 239 rfphon(2)=1 240 rfqvec(3)=1 241 242 write(msg, '(2a)' ) ch10," Extract the Phi^(1) coeficients from 3DTE" 243 call wrtout(std_out,msg,'COLL') 244 call ddb_lw%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec) 245 246 if (iblok == 0) then 247 call wrtout(std_out, " ") 248 call wrtout(std_out, "--- !WARNING") 249 call wrtout(std_out, sjoin("- Cannot find Phi^(1) tensor in DDB file:", filnamddb)) 250 call wrtout(std_out, " flexoflag=1 or 4 requires the DDB file to include the corresponding long wave 3rd derivatives") 251 end if 252 253 ! Look for the Gamma Block of the flexoelectric force response tensor in the DDB 254 qphon(:,:)=zero 255 qphnrm(:)=one 256 rfphon(:)=0 257 rfphon(1)=1 258 rfstrs(2)=3 259 rfqvec(3)=1 260 261 write(msg, '(2a)' ) ch10," Extract the FxE force response coeficients from 3DTE" 262 call wrtout(std_out,msg,'COLL') 263 call ddb_lw%get_block(jblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec) 264 265 if (jblok == 0) then 266 call wrtout(std_out, " ") 267 call wrtout(std_out, "--- !WARNING") 268 call wrtout(std_out, sjoin("- Cannot find FxE force response tensor in DDB file:", filnamddb)) 269 call wrtout(std_out, " flexoflag=1 or 4 requires the DDB file to include the corresponding long wave 3rd derivatives") 270 end if 271 272 ! Look for the stress tensor in the DDB 273 qphon(:,:)=zero 274 qphnrm(:)=one 275 rfphon(:)=0 276 rfstrs(:)=0 277 rfqvec(:)=0 278 rfstrs(4)=3 279 280 write(msg, '(2a)' ) ch10," Extract the stress tensor from 1DTE" 281 call wrtout(std_out,msg,'COLL') 282 call ddb%get_block(lblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,4) 283 284 if (lblok == 0) then 285 call wrtout(std_out, " ") 286 call wrtout(std_out, "--- !WARNING") 287 call wrtout(std_out, sjoin("- Cannot find stress tensor in DDB file:", filnamddb)) 288 call wrtout(std_out, " flexoflag=1 or 4 requires the DDB file to include the corresponding 2nd derivatives") 289 call wrtout(std_out, " to compute the Lagrange Elastic Tensor ") 290 end if 291 292 if (iblok/=0.and.jblok/=0) then 293 call dtlattflexo(ddb%amu,ddb%val(:,:,lblok),ddb_lw%val(:,:,jblok),ddb_lw%val(:,:,iblok),ddb_version,& 294 & intstrn,lattflexo,ddb%mpert,ddb%natom,crystal%ntypat,piezofr,prtvol,psinvdm,crystal%typat,crystal%ucvol,zeff) 295 end if 296 end if 297 298 !Merge the three contributions and print the total FxE tensor 299 totflexo(:,:,:,:)=ciflexo(:,:,:,:)+mixflexo(:,:,:,:)+lattflexo(:,:,:,:) 300 301 iwrite = ab_out > 0 302 if (iwrite) then 303 write(msg,'(3a)')ch10,' TOTAL flexoelectric tensor (units= nC/m) ',ch10 304 call wrtout([ab_out,std_out],msg,'COLL') 305 write(msg,*)' xx yy zz yz xz xy' 306 call wrtout([ab_out,std_out],msg,'COLL') 307 do ivar=1,6 308 elfd=alpha(ivar) 309 qvecd=beta(ivar) 310 write(msg,'(3x,a2,6f12.6)') voigt(ivar),totflexo(elfd,qvecd,1,1),totflexo(elfd,qvecd,2,2),totflexo(elfd,qvecd,3,3),& 311 totflexo(elfd,qvecd,2,3),totflexo(elfd,qvecd,1,3),totflexo(elfd,qvecd,1,2) 312 call wrtout([ab_out,std_out],msg,'COLL') 313 end do 314 do ivar=4,6 315 elfd=beta(ivar) 316 qvecd=alpha(ivar) 317 write(msg,'(3x,a2,6f12.6)') voigt(ivar+3),totflexo(elfd,qvecd,1,1),totflexo(elfd,qvecd,2,2),totflexo(elfd,qvecd,3,3),& 318 totflexo(elfd,qvecd,2,3),totflexo(elfd,qvecd,1,3),totflexo(elfd,qvecd,1,2) 319 call wrtout([ab_out,std_out],msg,'COLL') 320 end do 321 end if 322 323 324 DBG_EXIT("COLL") 325 326 end subroutine ddb_flexo
m_ddb/dm_psinv [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
dm_psinv
FUNCTION
Computes the pseudoinverse of the dynamical matrix (PRB 72, 035105 (2005)). This piece of code is copied from m_ddb_internalstr.F90.
INPUTS
asr= if /=0 acustic sume rule is imposed on the dynamical matrix blkval(2,3,mpert,3,mpert)= 2nd derivatives wrt two atom displacements (at least) d2asr(2,3,natom,3,natom)=ASR-correction iout=out file number mpert=maximum number of ipert natom=number of atoms in unit cell
OUTPUT
kmatrix(3*natom,3*natom) = array with the pseudo-inverse of dynamical matrix
SOURCE
1268 subroutine dm_psinv(asr,blkval,d2asr,iout,kmatrix,mpert,natom) 1269 1270 !Arguments ------------------------------- 1271 !scalars 1272 integer,intent(in) :: asr,iout,mpert,natom 1273 !arrays 1274 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert) 1275 real(dp),intent(in) :: d2asr(2,3,natom,3,natom) 1276 real(dp),intent(out) :: kmatrix(3*natom,3*natom) 1277 1278 !Local variables ------------------------- 1279 !scalars 1280 integer :: ii1,ipertA,ivarA 1281 integer :: ii2,ipertB,ivarB 1282 integer :: ier 1283 character(len=500) :: message 1284 !arrays 1285 real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom) 1286 real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2) 1287 real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3) 1288 real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom) 1289 real(dp) :: eigval(3*natom-3),eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3) 1290 real(dp) :: eigvecp(2,3*natom,3*natom) 1291 real(dp) :: zhpev1(2,2*3*natom-4) 1292 real(dp) :: zhpev1p(2,2*3*natom-1),zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2) 1293 real(dp) :: d2cart(2,3*natom,3*natom) 1294 1295 ! ********************************************************************* 1296 1297 DBG_ENTER("COLL") 1298 1299 d2cart = zero 1300 do ipertA=1,natom 1301 do ii1=1,3 1302 ivarA=ii1+3*(ipertA-1) 1303 do ipertB=1,natom 1304 do ii2=1,3 1305 ivarB=ii2+3*(ipertB-1) 1306 d2cart(1,ivarA,ivarB)=blkval(1,ii1,ipertA,ii2,ipertB) 1307 end do 1308 end do 1309 end do 1310 end do 1311 1312 !Eventually impose the acoustic sum rule 1313 !FIXME: this might depend on ifcflag: impose that it is 0 or generalize 1314 call asria_corr(asr,d2asr,d2cart,natom,natom) 1315 !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart) 1316 kmatrix = d2cart(1,:,:) 1317 Apmatr(:,:)=kmatrix(:,:) 1318 1319 !DEBUG 1320 !write(std_out,'(/,a,/)')'the force constant matrix' 1321 !do ivarA=1,3*natom 1322 !write(std_out,'(/)') 1323 !do ivarB=1,3*natom 1324 !write(std_out,'(es16.6)')kmatrix(ivarB,ivarA) 1325 !end do 1326 !end do 1327 !ENDDEBUG 1328 1329 Nmatr(:,:)=0.0_dp 1330 do ivarA=1,3*natom 1331 do ivarB=1,3*natom 1332 if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0)then 1333 Nmatr(ivarA,ivarB)=one 1334 end if 1335 if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1)then 1336 Nmatr(ivarA,ivarB)=one 1337 end if 1338 if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2)then 1339 Nmatr(ivarA,ivarB)=one 1340 end if 1341 end do 1342 end do 1343 1344 !DEBUG 1345 !do ivarA=1,3*natom 1346 !write(std_out,'(/)') 1347 !do ivarB=1,3*natom 1348 !write(std_out,'(es16.6)')Nmatr(ivarB,ivarA) 1349 !end do 1350 !end do 1351 !ENDDEBUG 1352 1353 !starting the pseudoinverting processes 1354 !then get the eigenvectors of the big matrix,give values to matrixBp 1355 Bpmatr=0.0_dp 1356 ii1=1 1357 do ivarA=1,3*natom 1358 do ivarB=1,ivarA 1359 Bpmatr(1,ii1)=Nmatr(ivarB,ivarA) 1360 ii1=ii1+1 1361 end do 1362 end do 1363 1364 !Bpmatr(2,:) is the imaginary part of the force matrix 1365 !then call the subroutines CHPEV and ZHPEV to get the eigenvectors 1366 call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier) 1367 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 1368 1369 !DEBUG 1370 !the eigenval and eigenvec 1371 !write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors' 1372 !do ivarA=1,3*natom 1373 !write(std_out,'(/)') 1374 !write(std_out,'(es16.6)')eigvalp(ivarA) 1375 !end do 1376 !do ivarA=1,3*natom 1377 !write(std_out,'(/)') 1378 !do ivarB=1,3*natom 1379 !write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA) 1380 !end do 1381 !end do 1382 !ENDDEBUG 1383 1384 !Then do the multiplication to get the reduced matrix,in two steps 1385 !After this the force constant matrix is decouple in two bloks, 1386 !acoustic and optical ones 1387 Cpmatr(:,:)=0.0_dp 1388 do ivarA=1,3*natom 1389 do ivarB=1,3*natom 1390 do ii1=1,3*natom 1391 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*Apmatr(ii1,ivarB) 1392 end do 1393 end do 1394 end do 1395 1396 Apmatr(:,:)=0.0_dp 1397 do ivarA=1,3*natom 1398 do ivarB=1,3*natom 1399 do ii1=1,3*natom 1400 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*eigvecp(1,ii1,ivarB) 1401 end do 1402 end do 1403 end do 1404 1405 !DEBUG 1406 !the blok diago 1407 !write(std_out,'(/,a,/)')'matrixAp' 1408 !do ivarA=1,3*natom 1409 !write(std_out,'(/)') 1410 !do ivarB=1,3*natom 1411 !write(std_out,'(es16.6)')Apmatr(ivarA,ivarB) 1412 !end do 1413 !end do 1414 !ENDDEBUG 1415 1416 !Check the last three eigenvalues whether too large or not 1417 ivarB=0 1418 do ivarA=3*natom-2,3*natom 1419 if (ABS(Apmatr(ivarA,ivarA))>tol6)then 1420 ivarB=1 1421 end if 1422 end do 1423 1424 if(ivarB==1)then 1425 write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,& 1426 & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,& 1427 & ' are too large at Gamma point.',ch10,& 1428 & ' Increase cutoff energy or k-points sampling.',ch10,& 1429 & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom) 1430 ABI_WARNING(message) 1431 call wrtout(iout,message,'COLL') 1432 end if 1433 1434 !Give the value of reduced matrix form Apmatr to Amatr 1435 do ivarA=1,3*natom-3 1436 do ivarB=1,3*natom-3 1437 Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB) 1438 end do 1439 end do 1440 1441 !Now the reduced matrix is in the matrixA, the convert it 1442 !first give the give the value of matixB from matrixA 1443 ii1=1 1444 do ivarA=1,3*natom-3 1445 do ivarB=1,ivarA 1446 Bmatr(1,ii1)=Amatr(ivarB,ivarA) 1447 ii1=ii1+1 1448 end do 1449 end do 1450 Bmatr(2,:)=0.0_dp 1451 1452 !Call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues 1453 call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier) 1454 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 1455 1456 !Check the unstable phonon modes, if the first is negative then print 1457 !warning message 1458 if(eigval(1)<-1.0*tol8)then 1459 write(message,'(9a)') ch10,& 1460 & ' --- !WARNING',ch10,& 1461 & ' Unstable eigenvalue detected in force constant matrix at Gamma point',ch10,& 1462 & ' The system under calculation is physically unstable.',ch10,& 1463 & ' ---',ch10 1464 call wrtout(std_out,message,'COLL') 1465 end if 1466 1467 !Do the matrix mutiplication to get pseudoinverse inverse matrix 1468 Cmatr(:,:)=0.0_dp 1469 Amatr(:,:)=0.0_dp 1470 do ivarA=1,3*natom-3 1471 Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA) 1472 end do 1473 1474 do ivarA=1,3*natom-3 1475 do ivarB=1,3*natom-3 1476 do ii1=1,3*natom-3 1477 Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*Cmatr(ii1,ivarB) 1478 end do 1479 end do 1480 end do 1481 1482 1483 !The second multiplication 1484 Cmatr(:,:)=0.0_dp 1485 do ivarA=1,3*natom-3 1486 do ivarB=1,3*natom-3 1487 do ii1=1,3*natom-3 1488 Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+ Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1) 1489 end do 1490 end do 1491 end do 1492 1493 !DEBUG 1494 !write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix' 1495 !do ivarA=1,3*natom 1496 !write(std_out,'(/)') 1497 !do ivarB=1,3*natom 1498 !write(std_out,'(es16.6)')Cmatr(ivarA,ivarB) 1499 !end do 1500 !end do 1501 !ENDDEBUG 1502 1503 !So now the inverse of the reduced matrix is in the matrixC 1504 !now do another mutilplication to get the pseudoinverse of the original 1505 Cpmatr(:,:)=0.0_dp 1506 Apmatr(:,:)=0.0_dp 1507 do ivarA=1,3*natom-3 1508 do ivarB=1,3*natom-3 1509 Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB) 1510 end do 1511 end do 1512 1513 !Now times the eigvecp 1514 do ivarA=1,3*natom 1515 do ivarB=1,3*natom 1516 do ii1=1,3*natom 1517 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*& 1518 & Cpmatr(ii1,ivarB) 1519 end do 1520 end do 1521 end do 1522 Cpmatr(:,:)=0.0_dp 1523 do ivarA=1,3*natom 1524 do ivarB=1,3*natom 1525 do ii1=1,3*natom 1526 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+ Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1) 1527 end do 1528 end do 1529 end do 1530 1531 !Now the inverse is in Cpmatr 1532 kmatrix(:,:)=Cpmatr(:,:) 1533 !transfer the inverse of k-matrix back to the k matrix 1534 !so now the inverse of k matrix is in the kmatrix 1535 !ending the part for pseudoinversing the K matrix 1536 1537 DBG_EXIT("COLL") 1538 1539 end subroutine dm_psinv
m_ddb/dtciflexo [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
dtciflexo
FUNCTION
Reads the Clamped Ion Flexoelectric 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 mpert =maximum number of ipert natom= number of atoms in unit cell ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files. ucvol= unit cell volume
OUTPUT
ciflexo(3,3,3,3) = type-II Clamped Ion Flexoelectric Tensor
SOURCE
350 subroutine dtciflexo(blkval,ddb_version,mpert,natom,ciflexo,ucvol) 351 352 !Arguments ------------------------------- 353 !scalars 354 integer,intent(in) :: ddb_version,mpert,natom 355 real(dp),intent(in) :: ucvol 356 !arrays 357 real(dp),intent(in) :: blkval(2,3*mpert*3*mpert*3*mpert) 358 real(dp),intent(out) :: ciflexo(3,3,3,3) 359 360 !Local variables ------------------------- 361 !scalars 362 integer,parameter :: cvrsio8=20100401 363 integer :: elfd,istrs,ivarA,ri,strsd,strsd1,strsd2,strst,qvecd 364 logical :: iwrite 365 real(dp) :: fac 366 character(len=500) :: msg 367 real(dp),parameter :: confac=e_Cb/Bohr_meter*1.d9 368 !arrays 369 integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/) 370 real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert) 371 character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/) 372 373 ! ********************************************************************* 374 375 DBG_ENTER("COLL") 376 377 d3cart(1,:,:,:,:,:,:) = reshape(blkval(1,:),shape = (/3,mpert,3,mpert,3,mpert/)) 378 d3cart(2,:,:,:,:,:,:) = reshape(blkval(2,:),shape = (/3,mpert,3,mpert,3,mpert/)) 379 380 !Define the factors to apply if DDB file has been created with the old version of 381 !the longwave driver. 382 if (ddb_version <= cvrsio8) then 383 fac=-two/ucvol 384 ri=2 385 else 386 fac=one/ucvol 387 ri=1 388 end if 389 390 !Extraction of the clamped-ion flexoelectric coeficients 391 do qvecd=1,3 392 do istrs=1,6 393 strsd1=alpha(istrs) 394 strsd2=beta(istrs) 395 strst=natom+3; if (istrs>3) strst=natom+4 396 strsd=istrs; if (istrs>3) strsd=istrs-3 397 do elfd=1,3 398 ciflexo(elfd,qvecd,strsd1,strsd2)=fac*d3cart(ri,elfd,natom+2,strsd,strst,qvecd,natom+8)*confac 399 if (istrs>3) ciflexo(elfd,qvecd,strsd2,strsd1)=ciflexo(elfd,qvecd,strsd1,strsd2) 400 end do 401 end do 402 end do 403 404 !Print results 405 iwrite = ab_out > 0 406 if (iwrite) then 407 write(msg,'(3a)')ch10,' Type-II electronic (clamped ion) flexoelectric tensor (units= nC/m) ',ch10 408 call wrtout([ab_out,std_out],msg,'COLL') 409 write(msg,*)' xx yy zz yz xz xy' 410 call wrtout([ab_out,std_out],msg,'COLL') 411 do ivarA=1,6 412 elfd=alpha(ivarA) 413 qvecd=beta(ivarA) 414 write(msg,'(3x,a2,6f12.6)') voigt(ivarA), ciflexo(elfd,qvecd,1,1),ciflexo(elfd,qvecd,2,2), & 415 ciflexo(elfd,qvecd,3,3), ciflexo(elfd,qvecd,2,3), & 416 ciflexo(elfd,qvecd,1,3),ciflexo(elfd,qvecd,1,2) 417 call wrtout([ab_out,std_out],msg,'COLL') 418 end do 419 do ivarA=4,6 420 elfd=beta(ivarA) 421 qvecd=alpha(ivarA) 422 write(msg,'(3x,a2,6f12.6)') voigt(ivarA+3), ciflexo(elfd,qvecd,1,1),ciflexo(elfd,qvecd,2,2), & 423 ciflexo(elfd,qvecd,3,3), ciflexo(elfd,qvecd,2,3), & 424 ciflexo(elfd,qvecd,1,3),ciflexo(elfd,qvecd,1,2) 425 call wrtout([ab_out,std_out],msg,'COLL') 426 end do 427 end if 428 429 DBG_EXIT("COLL") 430 431 end subroutine dtciflexo
m_ddb/dtlattflexo [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
dtlattflexo
FUNCTION
Reads Phi^(1), flexoelectric force response and internal strain tensors in the Gamma Block coming from the Derivative Data Base (long wave third-order derivatives). And computes the lattice mediated contribution to the flexoelectric tensor. It also computes and writes the Lagrangian Elastic tensor.
INPUTS
amu(ntypat)=mass each atom type in the unit cell blkval1d(2,3,mpert,3,mpert)= 1st derivative wrt stress (at least) blkval2d(2,3,mpert,3,mpert)= 2nd derivatives wrt atom displacements and electric field (at least) blkvalA(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies for FxE force response tensor blkvalB(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies for Phi^(1) tensor ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files. intstrn(3,3,3,natom)= relaxed-ion internal strain tensor mpert= maximum number of ipert natom= number of atoms in unit cell piezofr(3,natom,3,3)= piezoelectric force response tensor (required to compute the Lagrange elastic tensor) prtvol= if >1 print all tensors entering the structure of lattflexo psinvdm(3*natom,3*natom) = pseudo inverse of dynamical matrix typat(natom)= Type of each atom in the unit cell ucvol= unit cell volume
OUTPUT
lattflexo(3,3,3,3) = type-II lattice contribution to the Flexoelectric Tensor
SOURCE
706 subroutine dtlattflexo(amu,blkval1d,blkvalA,blkvalB,ddb_version,intstrn,lattflexo,mpert,natom,& 707 & ntypat,piezofr,prtvol,psinvdm,typat,ucvol,zeff) 708 709 !Arguments ------------------------------- 710 !scalars 711 integer,intent(in) :: ddb_version,mpert,natom,ntypat,prtvol 712 real(dp),intent(in) :: ucvol 713 714 !arrays 715 integer,intent(in) :: typat(natom) 716 real(dp),intent(in) :: amu(ntypat) 717 real(dp),intent(in) :: blkval1d(2,3,mpert,3,mpert) 718 real(dp),intent(in) :: blkvalA(2,3*mpert*3*mpert*3*mpert) 719 real(dp),intent(in) :: blkvalB(2,3*mpert*3*mpert*3*mpert) 720 real(dp),intent(in) :: intstrn(3,3,3,natom) 721 real(dp),intent(in) :: piezofr(3,natom,3,3) 722 real(dp),intent(in) :: psinvdm(3*natom,3*natom) 723 real(dp),intent(in) :: zeff(3,3,natom) 724 real(dp),intent(out) :: lattflexo(3,3,3,3) 725 726 !Local variables ------------------------- 727 !scalars 728 integer,parameter :: cvrsio8=20100401 729 integer :: elfd,iat,iatd,istrs,ivar,jat,jatd,jvar,kat,katd,strsd 730 integer :: strsd1,strsd2,strst,qvecd,qvecd2 731 logical :: iwrite 732 real(dp) :: fac,mtot 733 real(dp),parameter :: confac=e_Cb/Bohr_meter*1.d9 734 character(len=500) :: msg 735 !arrays 736 integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/) 737 real(dp) :: frcelast_t2(3,3,3,3) 738 real(dp) :: Csupkap(3,natom,3,3,3),Csupkapsum(3,3,3,3) 739 real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert) 740 real(dp) :: flexois(3,natom,3,3,3) 741 real(dp) :: flexofr(3,natom,3,3,3) 742 real(dp) :: hatCsupkap(3,natom,3,3,3) 743 real(dp) :: lmcelast(3,3,3,3) 744 real(dp) :: phi1(3,natom,3,natom,3) 745 real(dp) :: ricelast_t2(3,3,3,3) 746 real(dp) :: roundbkt_k(3,3,3,3,natom) 747 real(dp) :: sqrbkt_t1(3,3,3,3) 748 real(dp) :: stress(3,3) 749 character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/) 750 751 ! MR: Kept for testing 752 ! integer :: i,j,k,l 753 ! real(dp) :: delik,deljk,delil,deljl 754 755 ! ********************************************************************* 756 757 DBG_ENTER("COLL") 758 759 !Extraction of the stress tensor 760 stress(1,1)=blkval1d(1,1,natom+3,1,1) 761 stress(2,2)=blkval1d(1,2,natom+3,1,1) 762 stress(3,3)=blkval1d(1,3,natom+3,1,1) 763 stress(2,3)=blkval1d(1,1,natom+4,1,1);stress(3,2)=stress(2,3) 764 stress(1,3)=blkval1d(1,2,natom+4,1,1);stress(3,1)=stress(1,3) 765 stress(1,2)=blkval1d(1,3,natom+4,1,1);stress(2,1)=stress(1,2) 766 767 !Calculate the sublattice-dependent round bracket tensor of PRB 88,174106 (2013) 768 !First we need to extract the Phi^(1) tensor 769 d3cart(1,:,:,:,:,:,:) = reshape(blkvalB(1,:),shape = (/3,mpert,3,mpert,3,mpert/)) 770 d3cart(2,:,:,:,:,:,:) = reshape(blkvalB(2,:),shape = (/3,mpert,3,mpert,3,mpert/)) 771 772 !Define the factors to apply if DDB file has been created with the old version of 773 !the longwave driver. 774 if (ddb_version <= cvrsio8) then 775 fac=-two 776 else 777 fac=-one 778 end if 779 780 do qvecd=1,3 781 do jat=1,natom 782 do jatd=1,3 783 do iat=1,natom 784 do iatd=1,3 785 phi1(iatd,iat,jatd,jat,qvecd)=fac*d3cart(2,iatd,iat,jatd,jat,qvecd,natom+8) 786 end do 787 end do 788 end do 789 end do 790 end do 791 792 !Now perform the multiplication with the internal strain 793 roundbkt_k(:,:,:,:,:)=zero 794 do iat=1,natom 795 do iatd=1,3 796 do qvecd2=1,3 797 do jatd=1,3 798 do qvecd=1,3 799 do katd=1,3 800 do kat=1,natom 801 roundbkt_k(iatd,qvecd,jatd,qvecd2,iat)=roundbkt_k(iatd,qvecd,jatd,qvecd2,iat) & 802 + phi1(iatd,iat,katd,kat,qvecd)*intstrn(qvecd2,jatd,katd,kat) 803 end do 804 end do 805 ! write(100,'(5i3,1x,f12.6)') iat, iatd,qvecd,jatd,qvecd2,roundbkt_k(iatd,qvecd,jatd,qvecd2,iat) 806 end do 807 end do 808 end do 809 end do 810 end do 811 812 !Calculate now the Lagrange elastic tensors 813 !First we need to extract the clamped-ion flexoelectric force response tensor 814 d3cart(1,:,:,:,:,:,:) = reshape(blkvalA(1,:),shape = (/3,mpert,3,mpert,3,mpert/)) 815 d3cart(2,:,:,:,:,:,:) = reshape(blkvalA(2,:),shape = (/3,mpert,3,mpert,3,mpert/)) 816 817 !Define the factors to apply if DDB file has been created with the old version of 818 !the longwave driver. 819 if (ddb_version <= cvrsio8) then 820 fac=-two 821 else 822 fac=one 823 end if 824 825 do istrs=1,6 826 strsd1=alpha(istrs) 827 strsd2=beta(istrs) 828 strst=natom+3; if (istrs>3) strst=natom+4 829 strsd=istrs; if (istrs>3) strsd=istrs-3 830 do qvecd=1,3 831 do iat=1,natom 832 do iatd=1,3 833 flexofr(iatd,iat,qvecd,strsd1,strsd2)=fac*d3cart(1,iatd,iat,strsd,strst,qvecd,natom+8) 834 if (istrs>3) flexofr(iatd,iat,qvecd,strsd2,strsd1)=flexofr(iatd,iat,qvecd,strsd1,strsd2) 835 end do 836 end do 837 end do 838 end do 839 840 !Now compute the type-II frozen-ion elastic tensor (without stress corrected) 841 frcelast_t2(:,:,:,:)=zero 842 fac=one/ucvol 843 do strsd2=1,3 844 do strsd1=1,3 845 do qvecd=1,3 846 do iatd=1,3 847 do iat=1,natom 848 frcelast_t2(iatd,qvecd,strsd1,strsd2)=frcelast_t2(iatd,qvecd,strsd1,strsd2) + & 849 & flexofr(iatd,iat,qvecd,strsd1,strsd2)*fac 850 end do 851 ! write(100,'(4i3,1x,f12.6)') iatd,qvecd,strsd1,strsd2, frcelast_t2(iatd,qvecd,strsd1,strsd2) 852 end do 853 end do 854 end do 855 end do 856 857 !Now convert to type-I to obtain the square bracketed tensor of Born and Huang 858 do qvecd=1,3 859 do strsd2=1,3 860 do strsd1=1,3 861 do iatd=1,3 862 sqrbkt_t1(iatd,strsd1,strsd2,qvecd)=half*(frcelast_t2(iatd,qvecd,strsd1,strsd2) + & 863 & frcelast_t2(iatd,strsd2,strsd1,qvecd)) 864 end do 865 end do 866 end do 867 end do 868 869 !Now correct the stress in the square bracketed tesnor tensor 870 do qvecd=1,3 871 do strsd2=1,3 872 do strsd1=1,3 873 do iatd=1,3 874 if (iatd==strsd1) then 875 sqrbkt_t1(iatd,strsd1,strsd2,qvecd)=sqrbkt_t1(iatd,strsd1,strsd2,qvecd) - stress(strsd2,qvecd) 876 endif 877 end do 878 end do 879 end do 880 end do 881 882 !Now convert back to type-II in order to obtain the frozen ion Lagrange elastic tensor. 883 frcelast_t2(:,:,:,:)=zero 884 do iatd=1,3 885 do qvecd=1,3 886 do strsd1=1,3 887 do strsd2=1,3 888 frcelast_t2(iatd,qvecd,strsd1,strsd2)=sqrbkt_t1(iatd,strsd1,qvecd,strsd2) + & 889 & sqrbkt_t1(iatd,strsd2,strsd1,qvecd)-sqrbkt_t1(iatd,qvecd,strsd2,strsd1) 890 end do 891 end do 892 end do 893 end do 894 895 !MR: kept for testing only. If uncommented the resulting elastic tensors must agree with HWRV's ones 896 ! write(ab_out,*) 897 ! write(ab_out,*) "Stress tensor" 898 ! do i=1,3 899 ! write(ab_out,*) stress(i,:) 900 ! end do 901 ! do i=1,3 902 ! do j=1,3 903 ! do k=1,3 904 ! delik=0.d0; if(i==k) delik=1.d0 905 ! deljk=0.d0; if(j==k) deljk=1.d0 906 ! do l=1,3 907 ! delil=0.d0; if(i==l) delil=1.d0 908 ! deljl=0.d0; if(j==l) deljl=1.d0 909 910 ! frcelast_t2(i,j,k,l)=frcelast_t2(i,j,k,l) + & 911 ! & 0.5d0*( deljk*stress(i,l) + delik*stress(j,l) + delil*stress(j,k) & 912 ! & + deljl*stress(i,k) ) 913 914 ! end do 915 ! end do 916 ! end do 917 ! end do 918 919 !Now compute the contribution to the elastic tensor due to ion relaxations 920 !and sum with the clamped ion elastic tensor to obtain the relaxed ion one 921 lmcelast(:,:,:,:)=zero 922 do strsd2=1,3 923 do strsd1=1,3 924 do qvecd=1,3 925 do jatd=1,3 926 do iatd=1,3 927 do iat=1,natom 928 lmcelast(jatd,qvecd,strsd1,strsd2)=lmcelast(jatd,qvecd,strsd1,strsd2) - & 929 intstrn(jatd,qvecd,iatd,iat)*piezofr(iatd,iat,strsd1,strsd2)*fac 930 end do 931 end do 932 ricelast_t2(jatd,qvecd,strsd1,strsd2)=frcelast_t2(jatd,qvecd,strsd1,strsd2) + & 933 & lmcelast(jatd,qvecd,strsd1,strsd2) 934 end do 935 end do 936 end do 937 end do 938 939 !In last place compute the lattice contribution to the FxE tensor 940 !First obtain the C^{\kappa} tensor of Eq. 59 of PRB 88,174106 (2013) 941 !and its sublattice summation 942 Csupkapsum(:,:,:,:)=zero 943 do strsd2=1,3 944 do strsd1=1,3 945 do qvecd=1,3 946 do iatd=1,3 947 do iat=1,natom 948 Csupkap(iatd,iat,qvecd,strsd1,strsd2)=flexofr(iatd,iat,qvecd,strsd1,strsd2) + & 949 & roundbkt_k(iatd,qvecd,strsd1,strsd2,iat) 950 Csupkapsum(iatd,qvecd,strsd1,strsd2)=Csupkapsum(iatd,qvecd,strsd1,strsd2) + & 951 & Csupkap(iatd,iat,qvecd,strsd1,strsd2) 952 end do 953 end do 954 end do 955 end do 956 end do 957 958 !Then separate the mass-dependent part 959 mtot=zero 960 do iat=1,natom 961 mtot=mtot + amu(typat(iat)) 962 end do 963 964 do strsd2=1,3 965 do strsd1=1,3 966 do qvecd=1,3 967 do iatd=1,3 968 do iat=1,natom 969 hatCsupkap(iatd,iat,qvecd,strsd1,strsd2)=Csupkap(iatd,iat,qvecd,strsd1,strsd2) - & 970 & amu(typat(iat))/mtot*Csupkapsum(iatd,qvecd,strsd1,strsd2) 971 end do 972 end do 973 end do 974 end do 975 end do 976 977 !Now compute the type-II flexoelectric internal strain tensor 978 flexois(:,:,:,:,:)=zero 979 do strsd2=1,3 980 do strsd1=1,3 981 do qvecd=1,3 982 do iat=1,natom 983 do iatd=1,3 984 ivar=(iat-1)*3+iatd 985 do jat=1,natom 986 do jatd=1,3 987 jvar=(jat-1)*3+jatd 988 flexois(iatd,iat,qvecd,strsd1,strsd2)=flexois(iatd,iat,qvecd,strsd1,strsd2) + & 989 & psinvdm(ivar,jvar)*hatCsupkap(jatd,jat,qvecd,strsd1,strsd2) 990 end do 991 end do 992 end do 993 end do 994 end do 995 end do 996 end do 997 998 !Finally multiply by the effective charges to obtain the FxE tensor 999 lattflexo(:,:,:,:)=zero 1000 do strsd2=1,3 1001 do strsd1=1,3 1002 do qvecd=1,3 1003 do elfd=1,3 1004 do iat=1,natom 1005 do iatd=1,3 1006 lattflexo(elfd,qvecd,strsd1,strsd2)=lattflexo(elfd,qvecd,strsd1,strsd2) + & 1007 & zeff(elfd,iatd,iat)*flexois(iatd,iat,qvecd,strsd1,strsd2)/ucvol*confac 1008 end do 1009 end do 1010 end do 1011 end do 1012 end do 1013 end do 1014 1015 1016 iwrite = ab_out > 0 1017 if (iwrite) then 1018 write(msg,'(3a)')ch10,' Lagrange elastic tensor from long wave magnitudes (clamped ion) (units= 10^2 GPa) ',ch10 1019 call wrtout([ab_out,std_out],msg,'COLL') 1020 write(msg,*)' xx yy zz yz xz xy' 1021 call wrtout([ab_out,std_out],msg,'COLL') 1022 frcelast_t2(:,:,:,:)=frcelast_t2(:,:,:,:)*HaBohr3_GPa/100.00_dp 1023 do ivar=1,6 1024 iatd=alpha(ivar) 1025 qvecd=beta(ivar) 1026 write(msg,'(9f12.6)') frcelast_t2(iatd,qvecd,1,1), frcelast_t2(iatd,qvecd,2,2),frcelast_t2(iatd,qvecd,3,3), & 1027 frcelast_t2(iatd,qvecd,2,3), frcelast_t2(iatd,qvecd,1,3),frcelast_t2(iatd,qvecd,1,2) 1028 1029 call wrtout([ab_out,std_out],msg,'COLL') 1030 end do 1031 1032 write(msg,'(3a)')ch10,' Lagrange elastic tensor from long wave magnitudes (relaxed ion) (units= 10^2 GPa) ',ch10 1033 call wrtout([ab_out,std_out],msg,'COLL') 1034 write(msg,*)' xx yy zz yz xz xy' 1035 call wrtout([ab_out,std_out],msg,'COLL') 1036 ricelast_t2(:,:,:,:)=ricelast_t2(:,:,:,:)*HaBohr3_GPa/100.00_dp 1037 do ivar=1,6 1038 iatd=alpha(ivar) 1039 qvecd=beta(ivar) 1040 write(msg,'(9f12.6)') ricelast_t2(iatd,qvecd,1,1), ricelast_t2(iatd,qvecd,2,2),ricelast_t2(iatd,qvecd,3,3), & 1041 ricelast_t2(iatd,qvecd,2,3), ricelast_t2(iatd,qvecd,1,3),ricelast_t2(iatd,qvecd,1,2) 1042 1043 call wrtout([ab_out,std_out],msg,'COLL') 1044 end do 1045 if (prtvol > 1) then 1046 write(msg,'(3a)')ch10,' (...)^kappa contribution to the flexoelectric force-response tensor (units: eV)',ch10 1047 call wrtout([ab_out,std_out],msg,'COLL') 1048 write(msg,*)' atom dir xx yy zz yz xz xy' 1049 call wrtout([ab_out,std_out],msg,'COLL') 1050 roundbkt_k(:,:,:,:,:)=roundbkt_k(:,:,:,:,:)*Ha_eV 1051 do iat=1,natom 1052 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', roundbkt_k(1,1,1,1,iat), roundbkt_k(1,1,2,2,iat), & 1053 & roundbkt_k(1,1,3,3,iat),roundbkt_k(1,1,2,3,iat),roundbkt_k(1,1,1,3,iat),roundbkt_k(1,1,1,2,iat) 1054 call wrtout([ab_out,std_out],msg,'COLL') 1055 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', roundbkt_k(2,2,1,1,iat), roundbkt_k(2,2,2,2,iat), & 1056 & roundbkt_k(2,2,3,3,iat),roundbkt_k(2,2,2,3,iat),roundbkt_k(2,2,1,3,iat),roundbkt_k(2,2,1,2,iat) 1057 call wrtout([ab_out,std_out],msg,'COLL') 1058 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', roundbkt_k(3,3,1,1,iat), roundbkt_k(3,3,2,2,iat), & 1059 & roundbkt_k(3,3,3,3,iat),roundbkt_k(3,3,2,3,iat),roundbkt_k(3,3,1,3,iat),roundbkt_k(3,3,1,2,iat) 1060 call wrtout([ab_out,std_out],msg,'COLL') 1061 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', roundbkt_k(2,3,1,1,iat), roundbkt_k(2,3,2,2,iat), & 1062 & roundbkt_k(2,3,3,3,iat),roundbkt_k(2,3,2,3,iat),roundbkt_k(2,3,1,3,iat),roundbkt_k(2,3,1,2,iat) 1063 call wrtout([ab_out,std_out],msg,'COLL') 1064 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', roundbkt_k(1,3,1,1,iat), roundbkt_k(1,3,2,2,iat), & 1065 & roundbkt_k(1,3,3,3,iat),roundbkt_k(1,3,2,3,iat),roundbkt_k(1,3,1,3,iat),roundbkt_k(1,3,1,2,iat) 1066 call wrtout([ab_out,std_out],msg,'COLL') 1067 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', roundbkt_k(1,2,1,1,iat), roundbkt_k(1,2,2,2,iat), & 1068 & roundbkt_k(1,2,3,3,iat),roundbkt_k(1,2,2,3,iat),roundbkt_k(1,2,1,3,iat),roundbkt_k(1,2,1,2,iat) 1069 call wrtout([ab_out,std_out],msg,'COLL') 1070 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', roundbkt_k(3,2,1,1,iat), roundbkt_k(3,2,2,2,iat), & 1071 & roundbkt_k(3,2,3,3,iat),roundbkt_k(3,2,2,3,iat),roundbkt_k(3,2,1,3,iat),roundbkt_k(3,2,1,2,iat) 1072 call wrtout([ab_out,std_out],msg,'COLL') 1073 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', roundbkt_k(3,1,1,1,iat), roundbkt_k(3,1,2,2,iat), & 1074 & roundbkt_k(3,1,3,3,iat),roundbkt_k(3,1,2,3,iat),roundbkt_k(3,1,1,3,iat),roundbkt_k(3,1,1,2,iat) 1075 call wrtout([ab_out,std_out],msg,'COLL') 1076 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', roundbkt_k(2,1,1,1,iat), roundbkt_k(2,1,2,2,iat), & 1077 & roundbkt_k(2,1,3,3,iat),roundbkt_k(2,1,2,3,iat),roundbkt_k(2,1,1,3,iat),roundbkt_k(2,1,1,2,iat) 1078 call wrtout([ab_out,std_out],msg,'COLL') 1079 end do 1080 1081 write(msg,'(3a)')ch10,' [...]^kappa contribution to the flexoelectric force-response tensor (units: eV)',ch10 1082 call wrtout([ab_out,std_out],msg,'COLL') 1083 write(msg,*)' atom dir xx yy zz yz xz xy' 1084 call wrtout([ab_out,std_out],msg,'COLL') 1085 flexofr(:,:,:,:,:)=flexofr(:,:,:,:,:)*Ha_eV 1086 do iat=1,natom 1087 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', flexofr(1,iat,1,1,1),flexofr(1,iat,1,2,2),flexofr(1,iat,1,3,3),& 1088 & flexofr(1,iat,1,2,3),flexofr(1,iat,1,1,3),flexofr(1,iat,1,1,2) 1089 call wrtout([ab_out,std_out],msg,'COLL') 1090 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', flexofr(2,iat,2,1,1),flexofr(2,iat,2,2,2),flexofr(2,iat,2,3,3),& 1091 & flexofr(2,iat,2,2,3),flexofr(2,iat,2,1,3),flexofr(2,iat,2,1,2) 1092 call wrtout([ab_out,std_out],msg,'COLL') 1093 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', flexofr(3,iat,3,1,1),flexofr(3,iat,3,2,2),flexofr(3,iat,3,3,3),& 1094 & flexofr(3,iat,3,2,3),flexofr(3,iat,3,1,3),flexofr(3,iat,3,1,2) 1095 call wrtout([ab_out,std_out],msg,'COLL') 1096 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', flexofr(2,iat,3,1,1),flexofr(2,iat,3,2,2),flexofr(2,iat,3,3,3),& 1097 & flexofr(2,iat,3,2,3),flexofr(2,iat,3,1,3),flexofr(2,iat,3,1,2) 1098 call wrtout([ab_out,std_out],msg,'COLL') 1099 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', flexofr(1,iat,3,1,1),flexofr(1,iat,3,2,2),flexofr(1,iat,3,3,3),& 1100 & flexofr(1,iat,3,2,3),flexofr(1,iat,3,1,3),flexofr(1,iat,3,1,2) 1101 call wrtout([ab_out,std_out],msg,'COLL') 1102 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', flexofr(1,iat,2,1,1),flexofr(1,iat,2,2,2),flexofr(1,iat,2,3,3),& 1103 & flexofr(1,iat,2,2,3),flexofr(1,iat,2,1,3),flexofr(1,iat,2,1,2) 1104 call wrtout([ab_out,std_out],msg,'COLL') 1105 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', flexofr(3,iat,2,1,1),flexofr(3,iat,2,2,2),flexofr(3,iat,2,3,3),& 1106 & flexofr(3,iat,2,2,3),flexofr(3,iat,2,1,3),flexofr(3,iat,2,1,2) 1107 call wrtout([ab_out,std_out],msg,'COLL') 1108 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', flexofr(3,iat,1,1,1),flexofr(3,iat,1,2,2),flexofr(3,iat,1,3,3),& 1109 & flexofr(3,iat,1,2,3),flexofr(3,iat,1,1,3),flexofr(3,iat,1,1,2) 1110 call wrtout([ab_out,std_out],msg,'COLL') 1111 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', flexofr(2,iat,1,1,1),flexofr(2,iat,1,2,2),flexofr(2,iat,1,3,3),& 1112 & flexofr(2,iat,1,2,3),flexofr(2,iat,1,1,3),flexofr(2,iat,1,1,2) 1113 call wrtout([ab_out,std_out],msg,'COLL') 1114 end do 1115 end if 1116 1117 write(msg,'(3a)')ch10,' Flexoelectric force-response tensor (units: eV)',ch10 1118 call wrtout([ab_out,std_out],msg,'COLL') 1119 write(msg,*)' atom dir xx yy zz yz xz xy' 1120 call wrtout([ab_out,std_out],msg,'COLL') 1121 Csupkap(:,:,:,:,:)=Csupkap(:,:,:,:,:)*Ha_eV 1122 do iat=1,natom 1123 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', Csupkap(1,iat,1,1,1),Csupkap(1,iat,1,2,2),Csupkap(1,iat,1,3,3),& 1124 & Csupkap(1,iat,1,2,3),Csupkap(1,iat,1,1,3),Csupkap(1,iat,1,1,2) 1125 call wrtout([ab_out,std_out],msg,'COLL') 1126 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', Csupkap(2,iat,2,1,1),Csupkap(2,iat,2,2,2),Csupkap(2,iat,2,3,3),& 1127 & Csupkap(2,iat,2,2,3),Csupkap(2,iat,2,1,3),Csupkap(2,iat,2,1,2) 1128 call wrtout([ab_out,std_out],msg,'COLL') 1129 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', Csupkap(3,iat,3,1,1),Csupkap(3,iat,3,2,2),Csupkap(3,iat,3,3,3),& 1130 & Csupkap(3,iat,3,2,3),Csupkap(3,iat,3,1,3),Csupkap(3,iat,3,1,2) 1131 call wrtout([ab_out,std_out],msg,'COLL') 1132 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', Csupkap(2,iat,3,1,1),Csupkap(2,iat,3,2,2),Csupkap(2,iat,3,3,3),& 1133 & Csupkap(2,iat,3,2,3),Csupkap(2,iat,3,1,3),Csupkap(2,iat,3,1,2) 1134 call wrtout([ab_out,std_out],msg,'COLL') 1135 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', Csupkap(1,iat,3,1,1),Csupkap(1,iat,3,2,2),Csupkap(1,iat,3,3,3),& 1136 & Csupkap(1,iat,3,2,3),Csupkap(1,iat,3,1,3),Csupkap(1,iat,3,1,2) 1137 call wrtout([ab_out,std_out],msg,'COLL') 1138 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', Csupkap(1,iat,2,1,1),Csupkap(1,iat,2,2,2),Csupkap(1,iat,2,3,3),& 1139 & Csupkap(1,iat,2,2,3),Csupkap(1,iat,2,1,3),Csupkap(1,iat,2,1,2) 1140 call wrtout([ab_out,std_out],msg,'COLL') 1141 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', Csupkap(3,iat,2,1,1),Csupkap(3,iat,2,2,2),Csupkap(3,iat,2,3,3),& 1142 & Csupkap(3,iat,2,2,3),Csupkap(3,iat,2,1,3),Csupkap(3,iat,2,1,2) 1143 call wrtout([ab_out,std_out],msg,'COLL') 1144 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', Csupkap(3,iat,1,1,1),Csupkap(3,iat,1,2,2),Csupkap(3,iat,1,3,3),& 1145 & Csupkap(3,iat,1,2,3),Csupkap(3,iat,1,1,3),Csupkap(3,iat,1,1,2) 1146 call wrtout([ab_out,std_out],msg,'COLL') 1147 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', Csupkap(2,iat,1,1,1),Csupkap(2,iat,1,2,2),Csupkap(2,iat,1,3,3),& 1148 & Csupkap(2,iat,1,2,3),Csupkap(2,iat,1,1,3),Csupkap(2,iat,1,1,2) 1149 call wrtout([ab_out,std_out],msg,'COLL') 1150 end do 1151 if (prtvol > 1) then 1152 write(msg,'(3a)')ch10,' Flexoelectric force-response tensor minus mass dependent part (units: eV)',ch10 1153 call wrtout([ab_out,std_out],msg,'COLL') 1154 write(msg,*)' atom dir xx yy zz yz xz xy' 1155 call wrtout([ab_out,std_out],msg,'COLL') 1156 hatCsupkap(:,:,:,:,:)=hatCsupkap(:,:,:,:,:)*Ha_eV 1157 do iat=1,natom 1158 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'xx',hatCsupkap(1,iat,1,1,1),hatCsupkap(1,iat,1,2,2),hatCsupkap(1,iat,1,3,3),& 1159 & hatCsupkap(1,iat,1,2,3),hatCsupkap(1,iat,1,1,3),hatCsupkap(1,iat,1,1,2) 1160 call wrtout([ab_out,std_out],msg,'COLL') 1161 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'yy',hatCsupkap(2,iat,2,1,1),hatCsupkap(2,iat,2,2,2),hatCsupkap(2,iat,2,3,3),& 1162 & hatCsupkap(2,iat,2,2,3),hatCsupkap(2,iat,2,1,3),hatCsupkap(2,iat,2,1,2) 1163 call wrtout([ab_out,std_out],msg,'COLL') 1164 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'zz',hatCsupkap(3,iat,3,1,1),hatCsupkap(3,iat,3,2,2),hatCsupkap(3,iat,3,3,3),& 1165 & hatCsupkap(3,iat,3,2,3),hatCsupkap(3,iat,3,1,3),hatCsupkap(3,iat,3,1,2) 1166 call wrtout([ab_out,std_out],msg,'COLL') 1167 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'yz',hatCsupkap(2,iat,3,1,1),hatCsupkap(2,iat,3,2,2),hatCsupkap(2,iat,3,3,3),& 1168 & hatCsupkap(2,iat,3,2,3),hatCsupkap(2,iat,3,1,3),hatCsupkap(2,iat,3,1,2) 1169 call wrtout([ab_out,std_out],msg,'COLL') 1170 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'xz',hatCsupkap(1,iat,3,1,1),hatCsupkap(1,iat,3,2,2),hatCsupkap(1,iat,3,3,3),& 1171 & hatCsupkap(1,iat,3,2,3),hatCsupkap(1,iat,3,1,3),hatCsupkap(1,iat,3,1,2) 1172 call wrtout([ab_out,std_out],msg,'COLL') 1173 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'xy',hatCsupkap(1,iat,2,1,1),hatCsupkap(1,iat,2,2,2),hatCsupkap(1,iat,2,3,3),& 1174 & hatCsupkap(1,iat,2,2,3),hatCsupkap(1,iat,2,1,3),hatCsupkap(1,iat,2,1,2) 1175 call wrtout([ab_out,std_out],msg,'COLL') 1176 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'zy',hatCsupkap(3,iat,2,1,1),hatCsupkap(3,iat,2,2,2),hatCsupkap(3,iat,2,3,3),& 1177 & hatCsupkap(3,iat,2,2,3),hatCsupkap(3,iat,2,1,3),hatCsupkap(3,iat,2,1,2) 1178 call wrtout([ab_out,std_out],msg,'COLL') 1179 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'zx',hatCsupkap(3,iat,1,1,1),hatCsupkap(3,iat,1,2,2),hatCsupkap(3,iat,1,3,3),& 1180 & hatCsupkap(3,iat,1,2,3),hatCsupkap(3,iat,1,1,3),hatCsupkap(3,iat,1,1,2) 1181 call wrtout([ab_out,std_out],msg,'COLL') 1182 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'yx',hatCsupkap(2,iat,1,1,1),hatCsupkap(2,iat,1,2,2),hatCsupkap(2,iat,1,3,3),& 1183 & hatCsupkap(2,iat,1,2,3),hatCsupkap(2,iat,1,1,3),hatCsupkap(2,iat,1,1,2) 1184 call wrtout([ab_out,std_out],msg,'COLL') 1185 end do 1186 end if 1187 1188 write(msg,'(3a)')ch10,' Displacement-response flexoelectric internal strain tensor (units: Bohr^2)',ch10 1189 call wrtout([ab_out,std_out],msg,'COLL') 1190 write(msg,*)' atom dir xx yy zz yz xz xy' 1191 call wrtout([ab_out,std_out],msg,'COLL') 1192 do iat=1,natom 1193 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', flexois(1,iat,1,1,1),flexois(1,iat,1,2,2),flexois(1,iat,1,3,3),& 1194 & flexois(1,iat,1,2,3),flexois(1,iat,1,1,3),flexois(1,iat,1,1,2) 1195 call wrtout([ab_out,std_out],msg,'COLL') 1196 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', flexois(2,iat,2,1,1),flexois(2,iat,2,2,2),flexois(2,iat,2,3,3),& 1197 & flexois(2,iat,2,2,3),flexois(2,iat,2,1,3),flexois(2,iat,2,1,2) 1198 call wrtout([ab_out,std_out],msg,'COLL') 1199 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', flexois(3,iat,3,1,1),flexois(3,iat,3,2,2),flexois(3,iat,3,3,3),& 1200 & flexois(3,iat,3,2,3),flexois(3,iat,3,1,3),flexois(3,iat,3,1,2) 1201 call wrtout([ab_out,std_out],msg,'COLL') 1202 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', flexois(2,iat,3,1,1),flexois(2,iat,3,2,2),flexois(2,iat,3,3,3),& 1203 & flexois(2,iat,3,2,3),flexois(2,iat,3,1,3),flexois(2,iat,3,1,2) 1204 call wrtout([ab_out,std_out],msg,'COLL') 1205 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', flexois(1,iat,3,1,1),flexois(1,iat,3,2,2),flexois(1,iat,3,3,3),& 1206 & flexois(1,iat,3,2,3),flexois(1,iat,3,1,3),flexois(1,iat,3,1,2) 1207 call wrtout([ab_out,std_out],msg,'COLL') 1208 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', flexois(1,iat,2,1,1),flexois(1,iat,2,2,2),flexois(1,iat,2,3,3),& 1209 & flexois(1,iat,2,2,3),flexois(1,iat,2,1,3),flexois(1,iat,2,1,2) 1210 call wrtout([ab_out,std_out],msg,'COLL') 1211 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', flexois(3,iat,2,1,1),flexois(3,iat,2,2,2),flexois(3,iat,2,3,3),& 1212 & flexois(3,iat,2,2,3),flexois(3,iat,2,1,3),flexois(3,iat,2,1,2) 1213 call wrtout([ab_out,std_out],msg,'COLL') 1214 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', flexois(3,iat,1,1,1),flexois(3,iat,1,2,2),flexois(3,iat,1,3,3),& 1215 & flexois(3,iat,1,2,3),flexois(3,iat,1,1,3),flexois(3,iat,1,1,2) 1216 call wrtout([ab_out,std_out],msg,'COLL') 1217 write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', flexois(2,iat,1,1,1),flexois(2,iat,1,2,2),flexois(2,iat,1,3,3),& 1218 & flexois(2,iat,1,2,3),flexois(2,iat,1,1,3),flexois(2,iat,1,1,2) 1219 call wrtout([ab_out,std_out],msg,'COLL') 1220 end do 1221 1222 write(msg,'(3a)')ch10,' Type-II lattice contribution to flexoelectric tensor (units= nC/m) ',ch10 1223 call wrtout([ab_out,std_out],msg,'COLL') 1224 write(msg,*)' xx yy zz yz xz xy' 1225 call wrtout([ab_out,std_out],msg,'COLL') 1226 do ivar=1,6 1227 elfd=alpha(ivar) 1228 qvecd=beta(ivar) 1229 write(msg,'(3x,a2,6f12.6)') voigt(ivar),lattflexo(elfd,qvecd,1,1),lattflexo(elfd,qvecd,2,2),lattflexo(elfd,qvecd,3,3),& 1230 lattflexo(elfd,qvecd,2,3),lattflexo(elfd,qvecd,1,3),lattflexo(elfd,qvecd,1,2) 1231 call wrtout([ab_out,std_out],msg,'COLL') 1232 end do 1233 do ivar=4,6 1234 elfd=beta(ivar) 1235 qvecd=alpha(ivar) 1236 write(msg,'(3x,a2,6f12.6)') voigt(ivar+3),lattflexo(elfd,qvecd,1,1),lattflexo(elfd,qvecd,2,2),lattflexo(elfd,qvecd,3,3),& 1237 lattflexo(elfd,qvecd,2,3),lattflexo(elfd,qvecd,1,3),lattflexo(elfd,qvecd,1,2) 1238 call wrtout([ab_out,std_out],msg,'COLL') 1239 end do 1240 1241 end if 1242 DBG_EXIT("COLL") 1243 1244 end subroutine dtlattflexo
m_ddb/dtmixflexo [ Functions ]
[ Top ] [ m_ddb ] [ Functions ]
NAME
dtmixflexo
FUNCTION
Reads the P^(1) and Phi^(1) tensors in the Gamma Block coming from the Derivative Data Base (long wave third-order derivatives). And computes the mixed contribution to the flexoelectric tensor.
INPUTS
asr= if /=0 acustic sume rule is imposed on the dynamical matrix d2asr(2,3,natom,3,natom)=ASR-correction blkval1d(2,3,mpert,3,mpert)= 1st derivative wrt atom displacements (at least) blkval2d(2,3,mpert,3,mpert)= 2nd derivatives wrt two atom displacements (at least) blkval(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies for Phi^(1) tensor ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files. gprimd(3,3)= basis vectors in the reciprocal space intstrn_only= activates only the calculation of the internal strain tensor mpert =maximum number of ipert natom= number of atoms in unit cell piezofr(3,natom,3,3)= piezoelectric force response tensor pol1(3,3,3,natom)= tensor with the polarization induced by an atomic displacement (P^(1)) rprimd(3,3)= basis vectors in the real space ucvol= unit cell volume
OUTPUT
mixflexo(3,3,3,3) = type-II mixed contribution to the Flexoelectric Tensor intstrn(3,3,3,natom) = relaxed-ion internal strain tensor psinvdm(3*natom,3*natom) = pseudo inverse of dynamical matrix
SOURCE
467 subroutine dtmixflexo(asr,d2asr,blkval1d,blkval2d,blkval,ddb_version,gprimd,intstrn,intstrn_only, & 468 & mixflexo,mpert,natom,piezofr,pol1,psinvdm,rprimd,ucvol) 469 470 !Arguments ------------------------------- 471 !scalars 472 integer,intent(in) :: asr,ddb_version,mpert,natom 473 real(dp),intent(in) :: ucvol 474 logical,intent(in) :: intstrn_only 475 !arrays 476 real(dp),intent(in) :: d2asr(2,3,mpert,3,mpert) 477 real(dp),intent(in) :: blkval1d(2,3,mpert,3,mpert) 478 real(dp),intent(in) :: blkval2d(2,3,mpert,3,mpert) 479 real(dp),intent(in) :: blkval(2,3*mpert*3*mpert*3*mpert) 480 real(dp),intent(in) :: gprimd(3,3) 481 real(dp),intent(out) :: intstrn(3,3,3,natom) 482 real(dp),intent(out) :: piezofr(3,natom,3,3) 483 real(dp),intent(inout) :: pol1(3,3,3,natom) 484 real(dp),intent(out) :: psinvdm(3*natom,3*natom) 485 real(dp),intent(in) :: rprimd(3,3) 486 real(dp),intent(out) :: mixflexo(3,3,3,3) 487 488 !Local variables ------------------------- 489 !scalars 490 integer,parameter :: cvrsio8=20100401 491 integer :: elfd,iat,iatd,ivar,jat,jatd,jvar,katd,qvecd,qvecd2 492 logical :: iwrite 493 real(dp),parameter :: confac=e_Cb/Bohr_meter*1.d9 494 real(dp) :: fac 495 character(len=500) :: msg 496 !arrays 497 integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/) 498 real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert) 499 real(dp) :: redforces(3,natom),forces(3,natom) 500 real(dp) :: phi1(3,natom,3,natom,3) 501 integer :: flg1(3),flg2(3) 502 real(dp) :: vec1(3),vec2(3) 503 character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/) 504 505 ! ********************************************************************* 506 507 DBG_ENTER("COLL") 508 509 d3cart(1,:,:,:,:,:,:) = reshape(blkval(1,:),shape = (/3,mpert,3,mpert,3,mpert/)) 510 d3cart(2,:,:,:,:,:,:) = reshape(blkval(2,:),shape = (/3,mpert,3,mpert,3,mpert/)) 511 512 !P^(1) lacks the 1/ucvol factor 513 pol1=pol1/ucvol 514 515 !Define the factors to apply if DDB file has been created with the old version of 516 !the longwave driver. 517 if (ddb_version <= cvrsio8) then 518 fac=-two 519 else 520 fac=-one 521 end if 522 523 !Extraction of Phi^(1) tensor 524 do qvecd=1,3 525 do jat=1,natom 526 do jatd=1,3 527 do iat=1,natom 528 do iatd=1,3 529 phi1(iatd,iat,jatd,jat,qvecd)=fac*d3cart(2,iatd,iat,jatd,jat,qvecd,natom+8) 530 end do 531 end do 532 end do 533 end do 534 end do 535 536 !Extraction of the forces and conversion to cartesian coordinates 537 !to acount for the improper contribution 538 do iat=1,natom 539 do iatd=1,3 540 redforces(iatd,iat)=-blkval1d(1,iatd,iat,1,1) 541 end do 542 end do 543 forces(:,:)=redforces(:,:) 544 flg1(:)=1 545 do iat=1,natom 546 vec1(:)=forces(:,iat) 547 call cart39(flg1,flg2,gprimd,iat,natom,rprimd,vec1,vec2) 548 forces(:,iat)=vec2(:) 549 end do 550 551 !Calculate the piezoelectric force-response tensor including the improper contribution 552 piezofr(:,:,:,:)=zero 553 do qvecd=1,3 554 do iat=1,natom 555 do iatd=1,3 556 do jatd=1,3 557 do jat=1,natom 558 piezofr(iatd,iat,jatd,qvecd)=piezofr(iatd,iat,jatd,qvecd) + phi1(iatd,iat,jatd,jat,qvecd) 559 end do 560 if (iatd==qvecd) piezofr(iatd,iat,jatd,qvecd)=piezofr(iatd,iat,jatd,qvecd) + forces(jatd,iat) 561 end do 562 end do 563 end do 564 end do 565 566 !Calculate the ion-relaxed internal strain tensor 567 !First we need to obtain the pseudo-inverse of the dynamical matrix 568 call dm_psinv(asr,blkval2d,d2asr,ab_out,psinvdm,mpert,natom) 569 570 !Perfom the product with the piezo force-response 571 intstrn(:,:,:,:)=zero 572 do qvecd=1,3 573 do katd=1,3 574 do iatd=1,3 575 do iat=1,natom 576 ivar=(iat-1)*3+iatd 577 do jatd=1,3 578 do jat=1,natom 579 jvar=(jat-1)*3+jatd 580 581 intstrn(qvecd,katd,iatd,iat)= intstrn(qvecd,katd,iatd,iat) + & 582 psinvdm(ivar,jvar)*piezofr(jatd,jat,katd,qvecd) 583 584 end do 585 end do 586 end do 587 end do 588 end do 589 end do 590 591 if (.not.intstrn_only) then 592 !Finally calculate the mixed contribution to the FxE tensor 593 mixflexo(:,:,:,:)=zero 594 do elfd=1,3 595 do qvecd=1,3 596 do katd=1,3 597 do qvecd2=1,3 598 do iatd=1,3 599 do iat=1,natom 600 601 mixflexo(elfd,qvecd,katd,qvecd2)=mixflexo(elfd,qvecd,katd,qvecd2) - & 602 pol1(elfd,qvecd,iatd,iat)*intstrn(qvecd2,katd,iatd,iat)*confac 603 604 end do 605 end do 606 end do 607 end do 608 end do 609 end do 610 end if 611 612 !Print results 613 iwrite = ab_out > 0 614 if (iwrite) then 615 write(msg,'(3a)')ch10,' Force-response internal strain tensor from long-wave magnitudes (units: Hartree/Bohr)',ch10 616 call wrtout([ab_out,std_out],msg,'COLL') 617 write(msg,*)' atom dir xx yy zz yz xz xy' 618 call wrtout([ab_out,std_out],msg,'COLL') 619 do iat=1,natom 620 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'x', piezofr(1,iat,1,1),piezofr(1,iat,2,2),piezofr(1,iat,3,3),& 621 piezofr(1,iat,2,3),piezofr(1,iat,1,3),piezofr(1,iat,1,2) 622 call wrtout([ab_out,std_out],msg,'COLL') 623 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'y', piezofr(2,iat,1,1),piezofr(2,iat,2,2),piezofr(2,iat,3,3),& 624 piezofr(2,iat,2,3),piezofr(2,iat,1,3),piezofr(2,iat,1,2) 625 call wrtout([ab_out,std_out],msg,'COLL') 626 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'z', piezofr(3,iat,1,1),piezofr(3,iat,2,2),piezofr(3,iat,3,3),& 627 piezofr(3,iat,2,3),piezofr(3,iat,1,3),piezofr(3,iat,1,2) 628 call wrtout([ab_out,std_out],msg,'COLL') 629 end do 630 631 write(msg,'(3a)')ch10,' Displacement-response internal strain tensor from long-wave magnitudes (units: Bohr)',ch10 632 call wrtout([ab_out,std_out],msg,'COLL') 633 write(msg,*)' atom dir xx yy zz yz xz xy' 634 call wrtout([ab_out,std_out],msg,'COLL') 635 do iat=1,natom 636 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'x', intstrn(1,1,1,iat),intstrn(2,2,1,iat),intstrn(3,3,1,iat),& 637 intstrn(2,3,1,iat),intstrn(1,3,1,iat),intstrn(1,2,1,iat) 638 call wrtout([ab_out,std_out],msg,'COLL') 639 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'y', intstrn(1,1,2,iat),intstrn(2,2,2,iat),intstrn(3,3,2,iat),& 640 intstrn(2,3,2,iat),intstrn(1,3,2,iat),intstrn(1,2,2,iat) 641 call wrtout([ab_out,std_out],msg,'COLL') 642 write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'z', intstrn(1,1,3,iat),intstrn(2,2,3,iat),intstrn(3,3,3,iat),& 643 intstrn(2,3,3,iat),intstrn(1,3,3,iat),intstrn(1,2,3,iat) 644 call wrtout([ab_out,std_out],msg,'COLL') 645 end do 646 647 if (.not.intstrn_only) then 648 write(msg,'(3a)')ch10,' Type-II mixed contribution to flexoelectric tensor (units: nC/m)',ch10 649 call wrtout([ab_out,std_out],msg,'COLL') 650 write(msg,*)' xx yy zz yz xz xy' 651 call wrtout([ab_out,std_out],msg,'COLL') 652 do ivar=1,6 653 elfd=alpha(ivar) 654 qvecd=beta(ivar) 655 write(msg,'(3x,a2,6f12.6)') voigt(ivar),mixflexo(elfd,qvecd,1,1),mixflexo(elfd,qvecd,2,2),mixflexo(elfd,qvecd,3,3),& 656 mixflexo(elfd,qvecd,2,3),mixflexo(elfd,qvecd,1,3),mixflexo(elfd,qvecd,1,2) 657 call wrtout([ab_out,std_out],msg,'COLL') 658 end do 659 do ivar=4,6 660 elfd=beta(ivar) 661 qvecd=alpha(ivar) 662 write(msg,'(3x,a2,6f12.6)') voigt(ivar+3),mixflexo(elfd,qvecd,1,1),mixflexo(elfd,qvecd,2,2),mixflexo(elfd,qvecd,3,3),& 663 mixflexo(elfd,qvecd,2,3),mixflexo(elfd,qvecd,1,3),mixflexo(elfd,qvecd,1,2) 664 call wrtout([ab_out,std_out],msg,'COLL') 665 end do 666 end if 667 end if 668 669 DBG_EXIT("COLL") 670 671 end subroutine dtmixflexo