TABLE OF CONTENTS
ABINIT/update_e_field_vars [ Functions ]
NAME
update_e_field_vars
FUNCTION
This routine updates E field variables
COPYRIGHT
Copyright (C) 2003-2018 ABINIT group This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
atindx(natom)=index table for atoms, inverse of atindx (see gstate.f) atindx1(natom)=index table for atoms (see gstate.f) cg(2,mcg)=planewave coefficients of wavefunctions dimcprj(usepaw*natom)=lmn_size for each atom dtfil <type(datafiles_type)>=variables related to files gmet(3,3)=metric in reciprocal space gprimd(3,3)=reciprocal space dimensional primitive translations idir = determines directions for derivatives computed in ctocprj (0 for all) kg(3,mpw*mkmem)=reduced planewave coordinates mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mkmem=number of k points treated by this node. mpw=maximum dimensioned size of npw my_natom=number of atoms treated by current processor natom=number of atoms in cell nattyp(ntypat)=number of atoms of each type ngfft(18)=contain all needed information about 3D FFT, see ~ABINIT/Infos/vargs.htm#ngfft nkpt=number of k-points npwarr(nkpt)=number of planewaves in basis at this k point ntypat=number of types of atoms in unit cell pawrhoij(natom*usepaw) <type(pawrhoij_type)> atomic occupancies pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points (see initberry.f) pwind_alloc = first dimension of pwind pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations rmet(3,3)=metric in real space rprimd(3,3)=dimensional real space primitive translations (bohr) scfcv_level= 0 if calling before scf loop, 1 if during scfcv_quit=signals whether calling during scf quit (see scfcv.F90) scfcv_step=istep value of loop counter from scfcv.F90 ucvol=unit cell volume in bohr**3. unit_out= unit for output of the results (usually the .out file of ABINIT) The option unit_out = 0 is allowed. In this case, no information is written to the output file but only to the log file. usepaw= 1: use paw framework. 0:do not use paw. ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics
OUTPUT
efield_old_cart(3)=updating cartesian values of efield (used in berryopt 6,16,17) pel_cg(3)=electronic polarization pelev(3)=leading order PAW contribution in pel_cg (for reporting purposes only) pion(3)=ionic part of polarization ptot(3)=total polarization red_efield2=updating efield used in berryopt 16,17 red_efield2_old=updating efield used in berryopt 16.17 red_ptot=updating efield used in berryopt 16.17
SIDE EFFECTS
Input/Output dtset <type(dataset_type)>=all input variables in this dataset dtefield <type(efield_type)> = efield variables hdr <type(hdr_type)>=the header of wf, den and pot files mpi_enreg=information about MPI parallelization ptot_cart(3)=total polarization in cartesian coordinates xred(3,natom)=reduced atomic coordinates
TODO
NOTES
PARENTS
scfcv
CHILDREN
berryphase_new,ctocprj,getph,pawcprj_alloc,pawcprj_free,prtefield wrtout
SOURCE
92 #if defined HAVE_CONFIG_H 93 #include "config.h" 94 #endif 95 96 #include "abi_common.h" 97 98 subroutine update_e_field_vars(atindx,atindx1,cg,dimcprj,dtefield,dtfil,dtset,& 99 & efield_old_cart,gmet,gprimd,hdr,idir,kg,mcg,& 100 & mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,ngfft,nkpt,npwarr,ntypat,& 101 & pawrhoij,pawtab,pel_cg,pelev,pion,psps,ptot,ptot_cart,pwind,& 102 & pwind_alloc,pwnsfac,red_efield2,red_efield2_old,red_ptot,rmet,rprimd,& 103 & scfcv_level,scfcv_quit,scfcv_step,ucvol,unit_out,& 104 & usepaw,xred,ylm,ylmgr) 105 106 use defs_basis 107 use defs_datatypes 108 use defs_abitypes 109 use defs_wvltypes 110 use m_xmpi 111 use m_errors 112 use m_efield 113 use m_profiling_abi 114 115 use m_pawtab, only : pawtab_type 116 use m_pawrhoij, only : pawrhoij_type 117 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free 118 use m_kg, only : getph 119 120 !This section has been created automatically by the script Abilint (TD). 121 !Do not modify the following lines by hand. 122 #undef ABI_FUNC 123 #define ABI_FUNC 'update_e_field_vars' 124 use interfaces_14_hidewrite 125 use interfaces_66_nonlocal 126 use interfaces_67_common, except_this_one => update_e_field_vars 127 !End of the abilint section 128 129 implicit none 130 131 !Arguments ------------------------------------ 132 integer, intent(in) :: idir,mcg,mkmem,mpw,my_natom,natom,nkpt,ntypat 133 integer, intent(in) :: pwind_alloc,scfcv_level,scfcv_quit,scfcv_step,unit_out,usepaw 134 real(dp), intent(in) :: ucvol 135 type(datafiles_type), intent(in) :: dtfil 136 type(pseudopotential_type),intent(in) :: psps 137 type(dataset_type), intent(inout) :: dtset 138 type(efield_type), intent(inout) :: dtefield 139 type(hdr_type), intent(inout) :: hdr 140 type(MPI_type), intent(inout) :: mpi_enreg 141 !arrays 142 integer, intent(in) :: atindx(natom),atindx1(natom),dimcprj(usepaw*natom) 143 integer, intent(in) :: kg(3,mpw*mkmem),nattyp(ntypat) 144 integer, intent(in) :: ngfft(18),npwarr(nkpt),pwind(pwind_alloc,2,3) 145 real(dp), intent(in) :: cg(2,mcg),gmet(3,3),gprimd(3,3) 146 real(dp), intent(in) :: pwnsfac(2,pwind_alloc) 147 real(dp), intent(in) :: rmet(3,3),rprimd(3,3) 148 real(dp), intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm) 149 real(dp), intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm) 150 real(dp), intent(inout) :: ptot_cart(3),xred(3,natom),efield_old_cart(3) !vz_i 151 real(dp), intent(out) :: pel_cg(3),pelev(3),pion(3) !vz_i 152 real(dp), intent(inout) :: red_efield2(3),red_efield2_old(3) !vz_i 153 real(dp), intent(out) :: ptot(3),red_ptot(3) !vz_i 154 type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*usepaw) 155 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 156 157 !Local variables ------------------------- 158 !scalars 159 character(len=500) :: message 160 integer :: ctocprj_choice,iatom,ii,iorder_cprj,mcprj,my_nspinor,ncpgr 161 integer :: optberry,usecprj 162 logical :: calc_epaw3_force, calc_epaw3_stress, efield 163 !arrays 164 real(dp) :: efield_test_cart(3),red_efield1(3) 165 real(dp),allocatable :: ph1d(:,:) 166 type(pawcprj_type),allocatable :: cprj(:,:) 167 168 ! ************************************************************************* 169 170 efield = .false. 171 172 if ( dtset%berryopt == 4 .or. & 173 & dtset%berryopt == 6 .or. & 174 & dtset%berryopt == 7 .or. & 175 & dtset%berryopt ==14 .or. & 176 & dtset%berryopt ==16 .or. & 177 & dtset%berryopt ==17 ) efield = .true. 178 calc_epaw3_force = ( efield .and. dtset%optforces /= 0 .and. usepaw == 1 ) 179 calc_epaw3_stress = ( efield .and. dtset%optstress /= 0 .and. usepaw == 1 ) 180 181 usecprj=1; if (psps%usepaw==0) usecprj = 0 182 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 183 mcprj=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 184 185 ncpgr = 0 186 ctocprj_choice = 1 ! no derivs 187 if ( efield .and. psps%usepaw == 1) then 188 ABI_DATATYPE_ALLOCATE(cprj,(dtset%natom,mcprj)) 189 ! finite electric field may need gradients for forces, stress 190 if (calc_epaw3_force .and. .not. calc_epaw3_stress) then 191 ncpgr = 3; ctocprj_choice = 2 ! derivs w.r.t. position 192 else if (.not. calc_epaw3_force .and. calc_epaw3_stress) then 193 ncpgr = 6; ctocprj_choice = 3 ! derivs w.r.t strain 194 else if (calc_epaw3_force .and. calc_epaw3_stress) then 195 ncpgr = 9; ctocprj_choice = 23 ! derivs w.r.t. position and strain 196 end if 197 call pawcprj_alloc(cprj,ncpgr,dimcprj) 198 iatom=0 ; iorder_cprj=1 ! retain ordering of input list 199 ! all arguments to ctocprj are defined already except ph1d, do that here 200 ABI_ALLOCATE(ph1d,(2,3*(2*dtset%mgfft+1)*dtset%natom)) 201 call getph(atindx,dtset%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 202 call ctocprj(atindx,cg,ctocprj_choice,cprj,gmet,gprimd,iatom,idir,iorder_cprj,& 203 & dtset%istwfk,kg,dtset%kptns,mcg,mcprj,dtset%mgfft,dtset%mkmem,& 204 & mpi_enreg,psps%mpsang,dtset%mpw,dtset%natom,nattyp,dtset%nband,& 205 & dtset%natom,ngfft,dtset%nkpt,dtset%nloalg,npwarr,dtset%nspinor,& 206 & dtset%nsppol,dtset%ntypat,dtset%paral_kgb,ph1d,psps,rmet,& 207 & dtset%typat,ucvol,dtfil%unpaw,xred,ylm,ylmgr) 208 ABI_DEALLOCATE(ph1d) 209 else 210 ABI_DATATYPE_ALLOCATE(cprj,(0,0)) 211 end if ! end update of cprj 212 213 if ( efield ) then ! compute polarization and if necessary store cprj in efield 214 optberry=1 215 pel_cg(:) = zero;pelev=zero 216 call berryphase_new(atindx1,cg,cprj,dtefield,dtfil,dtset,psps,gprimd,hdr,psps%indlmn,kg,& 217 & psps%lmnmax,dtset%mband,mcg,mcprj,dtset%mkmem,mpi_enreg,dtset%mpw,my_natom,& 218 & dtset%natom,npwarr,dtset%nsppol,psps%ntypat,dtset%nkpt,optberry,pawrhoij,pawtab,& 219 & pel_cg,pelev,pion,ptot,red_ptot,pwind,& 220 & pwind_alloc,pwnsfac,rprimd,dtset%typat,ucvol,& 221 & unit_out,usecprj,psps%usepaw,xred,psps%ziontypat) 222 223 dtefield%red_ptot1(:)=red_ptot(:) 224 225 end if ! end compute polarization and store cprj for efield 226 227 if (efield .and. (scfcv_level == 0) ) then ! do this before scfcv loop 228 229 efield_old_cart(:)=dtset%efield(:) !!HONG 230 231 ! save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar 232 dtefield%efield2(:)=dtset%efield(:) 233 234 if ( dtset%berryopt ==16 .or. dtset%berryopt ==17) then !!HONG 235 do ii=1,3 236 red_efield2(ii)=zero 237 red_efield2_old(ii) =(ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) 238 end do 239 end if 240 241 if (dtset%berryopt == 14 .and. scfcv_quit /=1) then 242 ! ! Convert polarization to cartesian coords 243 244 ptot_cart(:)=zero 245 do ii = 1,3 246 ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + & 247 & rprimd(ii,3)*red_ptot(3) 248 end do 249 ptot_cart(:)=ptot_cart(:)/ucvol 250 251 do ii=1,3 252 dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) 253 end do 254 255 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 256 write(message,'(a,a)') ch10, 'scfcv: Constant reduced ebar-field:' 257 258 call wrtout(std_out,message,'COLL') 259 call prtefield(dtset,dtefield,std_out,rprimd) 260 261 if(dtset%prtvol>=10)then 262 call wrtout(ab_out,message,'COLL') 263 call prtefield(dtset,dtefield,ab_out,rprimd) 264 end if 265 266 ! updating E field 267 do ii =1,3 ! desired E field 268 efield_test_cart(ii)=gprimd(ii,1)*dtset%red_efieldbar(1) + & 269 & gprimd(ii,2)*dtset%red_efieldbar(2)+gprimd(ii,3)*dtset%red_efieldbar(3) 270 end do 271 272 ! if not convergence well, need to add some code here to make sure efield_test_cart(:) not change much 273 dtset%efield(:) = efield_test_cart(:) 274 275 end if ! berryopt ==14 276 277 end if ! end efield .and. scfcv_level 0 tasks 278 279 !!! 280 !!! Various printing and update steps for the different efield options 281 !!! 282 283 if (efield .and. (scfcv_level == 1) ) then ! do this each scf step 284 285 if (dtset%prtvol >= 10)then 286 write(message,'(6(a),3(e16.9,2x),a,a,3(e16.9,2x))')ch10,& 287 & ' scfcv: New value of the polarization:',ch10,& 288 & ' (reduced coordinates, a. u.)',ch10,& 289 & ' Electronic berry phase: ', (pel_cg(ii), ii = 1, 3) 290 call wrtout(ab_out,message,'COLL') 291 call wrtout(std_out,message,'COLL') 292 if(psps%usepaw==1) then 293 write(message,'(a,3(e16.9,2x))')& 294 & ' ...includes PAW on-site term: ', (pelev(ii), ii = 1, 3) 295 call wrtout(ab_out,message,'COLL') 296 call wrtout(std_out,message,'COLL') 297 end if 298 write(message,'(a,3(e16.9,2x),a,a,3(e16.9,2x))')& 299 & ' Ionic: ', (pion(ii), ii = 1, 3), ch10, & 300 & ' Total: ', (red_ptot(ii), ii = 1, 3) !!REC 301 call wrtout(ab_out,message,'COLL') 302 call wrtout(std_out,message,'COLL') 303 end if ! end prtvol >= 10 output 304 305 ptot_cart(:)=zero 306 do ii = 1,3 307 ptot_cart(ii)=rprimd(ii,1)*red_ptot(1) + rprimd(ii,2)*red_ptot(2) + & 308 & rprimd(ii,3)*red_ptot(3) 309 end do 310 ptot_cart(:)=ptot_cart(:)/ucvol 311 312 ! !=================================================================================================== 313 ! ! OUTPUT for fixed E 314 ! !=================================================================================================== 315 316 if (dtset%berryopt == 4) then 317 318 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 319 write(message,'(a,a)') ch10, 'scfcv: Constant unreduced E-field:' 320 call wrtout(std_out,message,'COLL') 321 call prtefield(dtset,dtefield,std_out,rprimd) 322 if(dtset%prtvol>=10)then 323 call wrtout(ab_out,message,'COLL') 324 call prtefield(dtset,dtefield,ab_out,rprimd) 325 end if 326 end if ! end berryopt 4 output 327 328 ! ===================================================================================== 329 ! ! fixed D calculation 330 ! !==================================================================================== 331 if (dtset%berryopt == 6) then 332 if (scfcv_step > 1) then 333 334 ! ! update efield taking damping into account dfield is in cartesian in dtset structure (contains input value) 335 ! ! same goes for efield - update the dtset%efield value 336 efield_test_cart(:)=dtset%ddamp*(dtset%dfield(:)-4.0d0*pi*ptot_cart(:))+& 337 & (1.0d0-dtset%ddamp)*efield_old_cart(:) 338 339 ! ! test whether change in efield in any direction exceed maxestep, if so, set the 340 ! ! change to maxestep instead ! need optimized ! 341 do ii = 1,3 342 343 if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then 344 345 write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii 346 write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), & 347 & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), & 348 & ", maxestep=",dtset%maxestep 349 350 351 if (efield_test_cart(ii) > efield_old_cart(ii)) then 352 efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep) 353 else 354 efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep) 355 end if 356 end if 357 end do 358 359 dtset%efield(:) = efield_test_cart(:) 360 361 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 362 write(message,'(a,a)') ch10, 'scfcv: Constant unreduced D-field - updating E-field:' 363 call wrtout(std_out,message,'COLL') 364 call prtefield(dtset,dtefield,std_out,rprimd) 365 if(dtset%prtvol>=10)then 366 call wrtout(ab_out,message,'COLL') 367 call prtefield(dtset,dtefield,ab_out,rprimd) 368 end if 369 370 ! ! need to update dtset%efield_dot(:) with new value 371 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 372 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 373 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 374 375 else 376 377 write(message,'(a,a)') ch10, 'scfcv: Constant unreduced D-field - Pre E-field:' 378 call wrtout(std_out,message,'COLL') 379 call prtefield(dtset,dtefield,std_out,rprimd) 380 if(dtset%prtvol>=10)then 381 call wrtout(ab_out,message,'COLL') 382 call prtefield(dtset,dtefield,ab_out,rprimd) 383 end if 384 385 end if ! scfcv_step >1 386 387 efield_old_cart(:)=dtset%efield(:) 388 end if ! berryopt ==6 389 ! !=================================================================================================== 390 ! ! fixed reduced d calculation 391 ! !=================================================================================================== 392 if (dtset%berryopt == 16) then 393 394 if (scfcv_step > 1) then 395 ! ! update efield taking damping into account reduced red_dfield 396 ! red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) 397 398 red_efield2(:)=dtset%ddamp*(dtset%red_dfield(:)-red_ptot(:))+ (1.0d0-dtset%ddamp)*red_efield2_old(:) 399 400 ! to calculate unreduced E 401 efield_test_cart(:)=(4*pi/ucvol)*(rprimd(:,1)*red_efield2(1)+rprimd(:,2)*red_efield2(2)+rprimd(:,3)*red_efield2(3)) 402 403 ! ! test whether change in efield in any direction exceed maxestep, if so, set the 404 ! ! change to maxestep instead ! need optimized ! 405 do ii = 1,3 406 407 if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then 408 409 write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii 410 write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), & 411 & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), & 412 & ", maxestep=",dtset%maxestep 413 414 if (efield_test_cart(ii) > efield_old_cart(ii)) then 415 efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep) 416 else 417 efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep) 418 end if 419 end if 420 end do 421 422 dtset%efield(:) = efield_test_cart(:) 423 424 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 425 write(message,'(a,a)') ch10, 'scfcv: Constant reduced d-field - updating E-field:' 426 call wrtout(std_out,message,'COLL') 427 call prtefield(dtset,dtefield,std_out,rprimd) 428 if(dtset%prtvol>=10)then 429 call wrtout(ab_out,message,'COLL') 430 call prtefield(dtset,dtefield,ab_out,rprimd) 431 end if 432 433 ! ! need to update dtset%efield_dot(:) with new value 434 ! ! This needs to be deleted when efield_dot is deleted 435 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 436 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 437 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 438 439 else 440 441 write(message,'(a,a)') ch10, 'scfcv: Constant reduced d-field - Pre E-field:' 442 call wrtout(std_out,message,'COLL') 443 call prtefield(dtset,dtefield,std_out,rprimd) 444 if(dtset%prtvol>=10)then 445 call wrtout(ab_out,message,'COLL') 446 call prtefield(dtset,dtefield,ab_out,rprimd) 447 end if 448 449 end if ! scfcv_step > 1 450 451 efield_old_cart(:)=dtset%efield(:) 452 red_efield2_old(:)=red_efield2(:) 453 end if ! berryopt ==16 454 455 456 ! !=================================================================================================== 457 ! ! fixed reduced d and ebar calculation (mixed BC) 458 ! !=================================================================================================== 459 if (dtset%berryopt == 17) then 460 461 if (scfcv_step > 1) then 462 ! ! update efield taking damping into account reduced red_dfield 463 ! red_efield1 and red_efield2 is reduced electric field, defined by Eq.(25) of Nat. Phys. suppl. (2009) 464 ! red_efield1 for fixed ebar, red_efield2 for fixed d calculation 465 466 ! save this value in order to print the final value of real electric field, comparing with the desired red_fieldbar 467 dtefield%efield2(:)=dtset%efield(:) 468 469 ! write(*,'(a,3i4)') "jfielddir=", (dtset%jfielddir(ii),ii=1,3) 470 471 do ii=1,3 472 if (dtset%jfielddir(ii) ==2 ) then ! direction under fixed d 473 dtset%red_efieldbar(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) ! update ebar which is not fixed 474 dtefield%efield_dot(ii) = dot_product(dtset%efield(:),rprimd(:,ii)) 475 red_efield2(ii)=dtset%ddamp*(dtset%red_dfield(ii) - red_ptot(ii)) + & 476 & (1.0d0-dtset%ddamp)*red_efield2_old(ii) ! d(ii) is fixed, update e(ii) may need ddamping here 477 478 ! write(message,'(a,a,i5,a,i5)') ch10, 'direction ', ii,' for fixed d, value is (2) ', dtset%jfielddir(ii) 479 ! call wrtout(ab_out,message,'COLL') 480 ! call wrtout(std_out,message,'COLL') 481 482 else if (dtset%jfielddir(ii) ==1 ) then ! direction under fixed ebar 483 red_efield2(ii)= (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) ! update e which is not fixed 484 dtset%red_dfield(ii)=red_ptot(ii) + (ucvol/(4*pi))*dot_product(dtset%efield(:),gprimd(:,ii)) ! update d 485 486 ! write(message,'(a,a,i5,a,i5)') ch10, 'direction ', ii,' for fixed ebar, value is (1) ', dtset%jfielddir(ii) 487 ! call wrtout(ab_out,message,'COLL') 488 ! call wrtout(std_out,message,'COLL') 489 490 end if 491 end do 492 493 do ii=1,3 494 red_efield1(ii) =(ucvol/(4*pi))*dot_product(dtset%red_efieldbar(:),gmet(:,ii)) 495 end do 496 497 498 dtset%red_efield(:)=(red_efield1(:) + red_efield2(:))/2.0d0 ! average reduced efield, 499 ! one is from fixed ebar part, 500 ! the other is from fixed d part. 501 ! This may need to be optimized !! 502 503 write(message,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Reduced efield from fixed ebar:', ch10, & 504 & ' e: ', (red_efield1(ii),ii=1,3), ch10 505 506 ! call wrtout(ab_out,message,'COLL') 507 call wrtout(std_out,message,'COLL') 508 509 write(message,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Reduced efield from fixed d:', ch10, & 510 & ' e: ', (red_efield2(ii),ii=1,3), ch10 511 512 ! call wrtout(ab_out,message,'COLL') 513 call wrtout(std_out,message,'COLL') 514 515 write(message,'(a,a,a,a,3(es16.9,2x),a)') ch10, 'Average reduced efield:', ch10, & 516 & ' e: ', (dtset%red_efield(ii),ii=1,3), ch10 517 518 ! call wrtout(ab_out,message,'COLL') 519 call wrtout(std_out,message,'COLL') 520 521 ! to calculate unreduced E 522 do ii=1,3 523 efield_test_cart(ii) = (4*pi/ucvol)* dot_product(dtset%red_efield(:),rprimd(:,ii)) 524 end do 525 526 ! ! test whether change in efield in any direction exceed maxestep, if so, set the 527 ! ! change to maxestep instead ! need optimized ! 528 do ii = 1,3 529 if (dabs(efield_test_cart(ii)-efield_old_cart(ii)) > dabs(dtset%maxestep)) then 530 531 write(std_out,'(a,a,i5)') "JH - "," E-field component:",ii 532 write(std_out,'(a,es13.5,a,es13.5,a,es13.5,a,es13.5)') " E(n)=",efield_test_cart(ii), & 533 & ", E(n-1)=",efield_old_cart(ii), ", E(n)-E(n-1)=", efield_test_cart(ii)-efield_old_cart(ii), & 534 & ", maxestep=",dtset%maxestep 535 536 if (efield_test_cart(ii) > efield_old_cart(ii)) then 537 efield_test_cart(ii) = efield_old_cart(ii) + dabs(dtset%maxestep) 538 else 539 efield_test_cart(ii) = efield_old_cart(ii) - dabs(dtset%maxestep) 540 end if 541 end if 542 end do 543 544 dtset%efield(:) = efield_test_cart(:) 545 546 ! !write the field parameters: D, E, P, d, e, p, dbar, ebar, pbar 547 write(message,'(a,a)') ch10, 'scfcv: Constant reduced ebar and d-field - updating E-field:' 548 call wrtout(std_out,message,'COLL') 549 call prtefield(dtset,dtefield,std_out,rprimd) 550 if(dtset%prtvol>=10)then 551 call wrtout(ab_out,message,'COLL') 552 call prtefield(dtset,dtefield,ab_out,rprimd) 553 end if 554 555 556 ! ! need to update dtset%efield_dot(:) with new value 557 ! ! This needs to be deleted when efield_dot is deleted 558 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 559 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 560 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 561 562 else 563 564 write(message,'(a,a)') ch10, 'scfcv: Constant reduced ebar and d-field - Pre E-field:' 565 call wrtout(std_out,message,'COLL') 566 call prtefield(dtset,dtefield,std_out,rprimd) 567 if(dtset%prtvol>=10)then 568 call wrtout(ab_out,message,'COLL') 569 call prtefield(dtset,dtefield,ab_out,rprimd) 570 end if 571 572 end if ! scfcv_step > 1 573 574 efield_old_cart(:)=dtset%efield(:) 575 red_efield2_old(:)=red_efield2(:) 576 577 end if ! berryopt ==17 578 579 end if ! end efield .and. scfcv_level 1 tasks 580 581 !deallocate cprj 582 if ( efield .and. psps%usepaw == 1) then 583 call pawcprj_free(cprj) 584 end if 585 ABI_DATATYPE_DEALLOCATE(cprj) 586 587 588 !DEBUG 589 !write(std_out,*)'update_e_field_vars exit' 590 !END_DEBUG 591 end subroutine update_e_field_vars