TABLE OF CONTENTS


ABINIT/update_e_field_vars [ Functions ]

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