TABLE OF CONTENTS


ABINIT/relaxpol [ Functions ]

[ Top ] [ Functions ]

NAME

 relaxpol

FUNCTION

 1) Compute polarization in cartesian coordinates
 2) Structural relaxation at fixed polarization: this routine
    solves the linear system of equations Eq.(13)
    of Na Sai et al., PRB 66, 104108 (2002).

COPYRIGHT

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

INPUTS

 blkflg(msize) = flag for every matrix element (0=> the element
   is not in the data block), (1=> the element is in the data blok)
 blkval(2,msize) = matrix that contains the second-order energy derivatives
 etotal = Kohn-Sham energy at zero electric field
 fred(3,natom) = -1 times the forces in reduced coordinates
 iatfix(natom) = indices of the atoms that are held fixed in the relaxation
 iout = unit number for output
 istrfix(6) = indices of the elements of the strain tensor that
   are held fixed in the relaxation
      1 = xx
      2 = yy
      3 = zz
      4 = yz & zy
      5 = xz & zx
      6 = xy & yx
 mpert = maximum number of ipert
 msize = dimension of blkflg and blkval
 natfix = number of atoms that are held fixed in the relaxation
 natom = number of atoms in the unit cell
 nstrfix = number of elements of the strain tensor that are held fixed in the relaxation
 pel(3) = electronic polarization not taking into account the factor 1/ucvol
  red_ptot(3) = total polarization reduced units   !!REC
 relaxat = 1: relax atomic positions
         = 0: do not relax atomic positions
 relaxstr = 1: relax cell parameters
          = 0: do not relax cell parameters
 strten(6) = stress tensor in cartesian coordinates
 targetpol(3) = target value of the polarization
 usepaw = 1 if PAW in use, 0 else (needed for polcart)

OUTPUT

NOTES

 - The elements of the dynamical matrix stored in blkval
   are symmetrized before computing the new atomic positions and cell parameters.
 - In case relaxat = 0 and relaxstr = 0, the routine only
   computes the polarization in cartesian coordinates.

PARENTS

      anaddb

CHILDREN

      dzgedi,dzgefa,matr3inv,polcart,symdyma,xcart2xred

SOURCE

 67 #if defined HAVE_CONFIG_H
 68 #include "config.h"
 69 #endif
 70 
 71 #include "abi_common.h"
 72 
 73 
 74 subroutine relaxpol(Crystal,blkflg,blkval,etotal,fred,iatfix,iout,istrfix,&
 75 & mpert,msize,natfix,natom,nstrfix,pel,red_ptot,relaxat,relaxstr,&
 76 & strten,targetpol,usepaw)
 77 
 78  use defs_basis
 79  use m_profiling_abi
 80  use m_errors
 81 
 82  use m_fstrings,  only : sjoin, itoa
 83  use m_abilasi,   only : dzgedi, dzgefa
 84  use m_geometry,  only : xcart2xred
 85  use m_dynmat,    only : symdyma
 86  use m_crystal,   only : crystal_t
 87 
 88 !This section has been created automatically by the script Abilint (TD).
 89 !Do not modify the following lines by hand.
 90 #undef ABI_FUNC
 91 #define ABI_FUNC 'relaxpol'
 92  use interfaces_32_util
 93  use interfaces_41_geometry
 94 !End of the abilint section
 95 
 96  implicit none
 97 
 98 !Arguments -------------------------------
 99 !scalars
100  integer,intent(in) :: iout,mpert,msize,natfix,natom,nstrfix
101  integer,intent(in) :: relaxat,relaxstr,usepaw
102  real(dp),intent(in) :: etotal
103  type(crystal_t),intent(in) :: Crystal
104 !arrays
105  integer,intent(in) :: blkflg(msize),iatfix(natom)
106  integer,intent(in) :: istrfix(6)
107  real(dp),intent(in) :: fred(3,natom),pel(3),strten(6)
108  real(dp),intent(in) :: red_ptot(3)
109  real(dp),intent(inout) :: blkval(2,msize),targetpol(3)
110 
111 !Local variables -------------------------
112 !scalars
113  integer :: flag,iatom,idir,ii,index,index1,index_tild,info,ipert,istrain
114  integer :: itypat,jdir,job,jpert,polunit,posi,posj,sizef
115  logical :: iwrite
116  real(dp) :: e1,fmax,poltmp,sigmax,tol,value,ucvol
117  character(len=500) :: message
118 !arrays
119  integer :: irelaxstrain(6)
120  integer,allocatable :: ipvt(:),irelaxat(:),rfpert(:,:)
121  real(dp) :: acell_new(3),delta_eta(6),delta_xcart(3,natom),det(2,2),diffpol(3),rprimd(3,3)
122  real(dp) :: diffsig(6),favg(3),gprimd(3,3),lambda(3),pel_cart(3),pelev(3)
123  real(dp) :: pion(3),pion_cart(3),ptot_cart(3),qphon(3),rprim(3,3)
124  real(dp) :: rprimd_new(3,3),sigelfd(6),strainmat(3,3)
125  real(dp) :: xcart_new(3,natom),xred_new(3,natom)
126  real(dp),allocatable :: cfac(:,:),delta(:),dymati(:),fcart(:,:),fcmat(:,:,:)
127  real(dp),allocatable :: fdiff(:,:),felfd(:,:),ifcmat(:,:,:),vec(:),zgwork(:,:)
128 
129 ! *********************************************************************
130 
131  rprimd = Crystal%rprimd
132  ucvol = Crystal%ucvol
133  iwrite = iout > 0
134 
135 !Check if some degrees of freedom remain fixed during the optimization
136 
137  ABI_ALLOCATE(irelaxat,(natom))
138  irelaxat(:) = 1   ; irelaxstrain(:) = 1
139  if (natfix > 0) then
140    do ii = 1, natfix
141      iatom = iatfix(ii)
142      if ((iatom > natom).or.(iatom < 0)) then
143        write(message, '(a,i0,a,i0,a,a,a,a,a)')&
144 &       'The value of iatfix(',ii,') is ',iatom,', which is not allowed.',ch10,&
145 &       'iatfix must be larger than 0 and smaller than natom.',ch10,&
146 &       'Action: correct iatfix in your input file.'
147        MSG_ERROR(message)
148      end if
149      irelaxat(iatom) = 0
150    end do
151  end if
152 
153  if (nstrfix > 0) then
154    do ii = 1, nstrfix
155      istrain = istrfix(ii)
156      if ((istrain > 6).or.(istrain < 0)) then
157        write(message, '(a,i0,a,i0,a,a,a,a,a)')&
158 &       'istrfix(',ii,') is',istrain,', which is not allowed.',ch10,&
159 &       'istrfix must be larger than 0 and smaller than 6.',ch10,&
160 &       'Action : correct istrfix in your input file.'
161        MSG_ERROR(message)
162      end if
163      irelaxstrain(istrain) = 0
164    end do
165  end if
166 
167 
168  ABI_ALLOCATE(rfpert,(mpert,3))
169  ABI_ALLOCATE(cfac,(mpert,mpert))
170  call matr3inv(rprimd,gprimd)
171 
172 !Compute the size of the matrix that contains the second-order derivatives
173 
174  sizef = 3
175  rfpert(:,:) = 0
176  rfpert(natom+2,1:3) = 1
177  if (relaxat == 1) then
178    do iatom = 1, natom
179      if (irelaxat(iatom) == 1) then
180        sizef = sizef + 3
181        rfpert(iatom,1:3) = 1
182      end if
183    end do
184  end if
185  ii = natom + 2
186  if (relaxstr == 1) then
187    istrain = 0
188    do ipert = (natom+3), (natom+4)
189      do idir = 1, 3
190        istrain = istrain + 1
191        if (irelaxstrain(istrain) == 1) then
192          sizef = sizef + 1
193          rfpert(ipert,idir) = 1
194        end if
195      end do
196    end do
197  end if
198 
199  ABI_ALLOCATE(fcmat,(2,sizef,sizef))
200  ABI_ALLOCATE(ifcmat,(2,sizef,sizef))
201  ABI_ALLOCATE(vec,(sizef))
202  ABI_ALLOCATE(delta,(sizef))
203  ABI_ALLOCATE(ipvt,(sizef))
204  ABI_ALLOCATE(zgwork,(2,sizef))
205  ABI_ALLOCATE(fcart,(3,natom))
206  ABI_ALLOCATE(felfd,(3,natom))
207  ABI_ALLOCATE(fdiff,(3,natom))
208 
209 !Build the vector that stores the forces, sigma and the polarization
210 
211  vec(:) = zero
212  posi = 0
213 
214  if (relaxat == 1) then
215 
216 !  Note conversion to cartesian coordinates (bohr) AND
217 !  negation to make a force out of a gradient
218 !  Also subtract off average force from each force
219 !  component to avoid spurious drifting of atoms across cell.
220    favg(:) = zero
221    do iatom = 1, natom
222      do idir = 1, 3
223        fcart(idir,iatom) = -(gprimd(idir,1)*fred(1,iatom) + &
224 &       gprimd(idir,2)*fred(2,iatom) + &
225 &       gprimd(idir,3)*fred(3,iatom))
226        favg(idir) = favg(idir) + fcart(idir,iatom)
227      end do
228    end do
229    favg(:) = favg(:)/dble(natom)
230    do iatom = 1, natom
231      fcart(:,iatom) = fcart(:,iatom) - favg(:)
232    end do
233 
234    do iatom = 1, natom
235      if (irelaxat(iatom) == 1) then
236        do idir = 1, 3
237          posi = posi + 1
238          vec(posi) = fcart(idir,iatom)
239        end do
240      end if
241    end do
242 
243  end if    ! relaxat == 1
244 
245 !DEBUG
246 !write(std_out,*)'Forces in cartesian coords'
247 !do iatom = 1, natom
248 !write(std_out,'(3(2x,e16.9))')(fcart(idir,iatom),idir = 1, 3)
249 !end do
250 !stop
251 !ENDDEBUG
252 
253 !Transform target polarization to atomic units
254  targetpol(:) = targetpol(:)*((Bohr_Ang*1.0d-10)**2)/e_Cb
255 
256 !Compute ionic polarization
257  pion(:) = zero
258  do iatom = 1, natom
259    itypat = Crystal%typat(iatom)
260    do idir = 1, 3
261      poltmp = Crystal%zion(itypat) * Crystal%xred(idir,iatom)
262      poltmp = poltmp - two*nint(poltmp/two)   ! fold into [-1,1]
263      pion(idir) = pion(idir) + poltmp
264    end do
265  end do
266  do idir = 1, 3
267    pion(idir) = pion(idir) - two*nint(pion(idir)/two) ! fold into [-1,1]
268  end do
269 
270 !Transform the polarization to cartesian coordinates
271  polunit = 3
272  pelev=zero
273  call polcart(red_ptot,pel,pel_cart,pelev,pion,pion_cart,polunit,&
274 & ptot_cart,rprimd,ucvol,iout,usepaw)
275 
276  do idir = 1, 3
277    posi = posi + 1
278    vec(posi) = ptot_cart(idir) - targetpol(idir)
279  end do
280 
281 
282  if (relaxstr == 1) then
283    do istrain = 1, 6
284      if (irelaxstrain(istrain) == 1) then
285        posi = posi + 1
286        vec(posi) = -1._dp*strten(istrain)*ucvol
287      end if
288    end do
289  end if
290 
291 
292 !Symmetrize the dynamical matrix
293 
294  ABI_ALLOCATE(dymati,(2*3*natom*3*natom))
295 !by the symdyma routine
296  do ipert = 1, natom
297    do idir = 1, 3
298      do jpert = 1, natom
299        do jdir = 1, 3
300          index  = jdir +3*((jpert - 1) + mpert*((idir - 1) + 3*(ipert - 1)))
301          index1 = jdir +3*((jpert - 1) + natom*((idir - 1) + 3*(ipert - 1)))
302          dymati(2*index1 - 1) = blkval(1,index)
303          dymati(2*index1    ) = blkval(2,index)
304        end do
305      end do
306    end do
307  end do
308 
309  qphon(:) = zero
310  call symdyma(dymati,Crystal%indsym,natom,Crystal%nsym,qphon,rprimd,Crystal%symrel,Crystal%symafm)
311 
312  do ipert = 1, natom
313    do idir = 1, 3
314      do jpert = 1, natom
315        do jdir = 1, 3
316          index  = jdir +3*((jpert - 1) + mpert*((idir - 1) + 3*(ipert - 1)))
317          index1 = jdir +3*((jpert - 1) + natom*((idir - 1) + 3*(ipert - 1)))
318          blkval(1,index) = dymati(2*index1 - 1)
319          blkval(2,index) = dymati(2*index1    )
320        end do
321      end do
322    end do
323  end do
324 
325  ABI_DEALLOCATE(dymati)
326 
327 !Define conversion factors for blkval
328  cfac(:,:) = 1._dp
329  cfac(1:natom,natom+2) = -1._dp/ucvol
330  cfac(natom+2,1:natom) = -1._dp/ucvol
331  cfac(natom+3:natom+4,natom+2) = -1._dp
332  cfac(natom+2,natom+3:natom+4) = -1._dp
333 
334 
335 !Build the matrix that contains the second-order derivatives
336 !ipert = natom + 1 corresponds to the ddk perturbation, that
337 !is not needed; so skip it
338 
339  fcmat(:,:,:) = zero
340 
341  posi = 0
342  flag = 0
343 ! When fcmat has been build, flag = 0 if all elements were available.
344 ! Otherwise, it will be 1. In case one element is missing, check if
345 ! it can be obtained by changing the order of the perturbations
346 
347  do ipert = 1, mpert
348    do idir = 1, 3
349      if (rfpert(ipert,idir) == 1) then
350        posi = posi + 1
351        posj = 0
352 
353        do jpert = 1, mpert
354          do jdir = 1, 3
355            if (rfpert(jpert,jdir) == 1) then
356              index = jdir +3*((jpert - 1) + mpert*((idir - 1) + 3*(ipert - 1)))
357              index_tild = idir +3*((ipert - 1) + mpert*((jdir - 1) + 3*(jpert - 1)))
358              posj = posj + 1
359              if ((ipert /= natom + 2).or.(jpert /= natom + 2)) then
360                if (blkflg(index) == 1) then
361                  fcmat(:,posi,posj) = blkval(:,index)*cfac(ipert,jpert)
362                else if (blkflg(index_tild) == 1) then
363                  fcmat(:,posi,posj) = blkval(:,index_tild)*cfac(ipert,jpert)
364                  blkval(:,index) = blkval(:,index_tild)
365                else
366                  flag = 1
367                  write(std_out,'(a,4(2x,i3))')'relaxpol: could not find element:',idir,ipert,jdir,jpert
368                end if
369              end if
370 !            DEBUG
371 !            write(100,'(4(2x,i3),5x,f16.9)')idir,ipert,jdir,jpert,fcmat(1,posi,posj)
372 !            ENDDEBUG
373            end if
374          end do
375        end do
376 
377      end if
378    end do
379  end do
380 
381  if (flag == 1) then
382    write(message, '(a,a,a,i0,a,i0,a,a,a,a)' )&
383 &   'Some of the second order derivatives required to deal with the case',ch10,&
384 &   'relaxat = ',relaxat,', relaxstr = ', relaxstr, ch10,&
385 &   'are missing in the DDB.',ch10,&
386 &   'Action: correct your DDB or change your input file.'
387    MSG_ERROR(message)
388  end if
389 
390 
391 !Compute the inverse of the force constant matrix
392 
393  if ((relaxat /= 0).or.(relaxstr /= 0)) then
394 
395    job = 1          ! compute inverse only
396    ifcmat(:,:,:) = fcmat(:,:,:)
397 
398    call dzgefa(ifcmat,sizef,sizef,ipvt,info)
399    ABI_CHECK(info == 0, sjoin("dzgefa returned:", itoa(info)))
400    call dzgedi(ifcmat,sizef,sizef,ipvt,det,zgwork,job)
401 
402 !  DEBUG
403 !  write(100,*)'relaxat = ',relaxat
404 !  write(100,*)'relaxstr = ',relaxstr
405 !  write(100,*)'irelaxat = '
406 !  write(100,*)irelaxat(:)
407 !  write(100,*)'irelaxstrain = '
408 !  write(100,*)irelaxstrain(:)
409 !  write(100,*)'sizef = ',sizef
410 !  write(100,*)'targetpol ='
411 !  write(100,*)targetpol(:)
412 !  do ipert = 1, sizef
413 !  do jpert = 1, sizef
414 !  write(100,'(2(2x,i3),2x,e16.9)')ipert,jpert,fcmat(1,ipert,jpert)
415 !  end do
416 !  end do
417 !  stop
418 !  ENDDEBUG
419 
420 !  Compute \delta R, \delta \eta and \lambda
421    delta(:) = zero
422    do ipert = 1, sizef
423      do jpert = 1, sizef
424        delta(ipert) = delta(ipert) + ifcmat(1,ipert,jpert)*vec(jpert)
425      end do
426    end do
427 
428 
429 !  Update atomic positions
430    posi = 0
431    if (relaxat == 1) then
432 
433      delta_xcart(:,:) = zero
434      xcart_new(:,:) = zero
435      do iatom = 1, natom
436        if (irelaxat(iatom) == 1) then
437          do idir = 1, 3
438            posi = posi + 1
439            delta_xcart(idir,iatom) = delta(posi)
440          end do
441        end if
442      end do
443 
444 !    Drop unsignificant digits in order to eleminate numerical noise
445      tol = 10000000._dp
446      do iatom = 1, natom
447        do idir = 1, 3
448          value = delta_xcart(idir,iatom)
449          ii = log10(abs(value))
450          if (ii <= 0) then
451            ii = abs(ii) + 1
452            value = one*int(tol*value*10.0_dp**ii)/(tol*10.0_dp**ii) !vz_d
453          else
454            value = one*int(tol*value/(10.0_dp**ii))*(10.0_dp**ii)/tol !vz_d
455          end if
456          delta_xcart(idir,iatom) = value
457        end do
458      end do
459 
460      xcart_new(:,:) = Crystal%xcart(:,:) + delta_xcart(:,:)
461      call xcart2xred(natom,rprimd,xcart_new,xred_new)
462    end if  ! relaxat == 1
463 
464 ! Compute lambda and the value of the energy functional F - \lambda \cdot P$
465 
466    e1 = etotal
467    do idir = 1, 3
468      posi = posi + 1
469      lambda(idir) = delta(posi)
470      e1 = e1 - lambda(idir)*ptot_cart(idir)
471    end do
472 
473 !  Update cell parameters
474    if (relaxstr == 1) then
475      delta_eta(:) = zero
476      do istrain = 1, 6
477        if (irelaxstrain(istrain) == 1) then
478          posi = posi + 1
479          delta_eta(istrain) = delta(posi)
480        end if
481      end do
482 
483      do istrain = 1, 3
484        strainmat(istrain,istrain) = delta_eta(istrain)
485      end do
486      strainmat(2,3) = delta_eta(4)/2._dp ; strainmat(3,2) = delta_eta(4)/2._dp
487      strainmat(1,3) = delta_eta(5)/2._dp ; strainmat(3,1) = delta_eta(5)/2._dp
488      strainmat(2,1) = delta_eta(6)/2._dp ; strainmat(1,2) = delta_eta(6)/2._dp
489 
490      rprimd_new(:,:) = 0._dp
491      do idir = 1, 3
492        do jdir = 1, 3
493          do ii = 1, 3
494            rprimd_new(jdir,idir) = rprimd_new(jdir,idir) + &
495 &           rprimd(ii,idir)*strainmat(ii,jdir)
496          end do
497        end do
498      end do
499      rprimd_new(:,:) = rprimd_new(:,:) + rprimd(:,:)
500 
501      acell_new(:) = zero
502      do idir = 1, 3
503        do jdir = 1, 3
504          acell_new(idir) = acell_new(idir) + &
505 &         rprimd_new(jdir,idir)*rprimd_new(jdir,idir)
506        end do
507        acell_new(idir) = sqrt(acell_new(idir))
508        rprim(:,idir) = rprimd_new(:,idir)/acell_new(idir)
509      end do
510 
511    end if          ! relaxstr == 1
512 
513 !  Write out the results
514 
515    if (iwrite) then
516      write(iout,*)
517      write(iout,'(a,80a,a)') ch10,('=',ii=1,80),ch10
518      write(iout,*)
519      write(iout,*)'Relaxation of the geometry at fixed polarization:'
520      write(iout,*)
521      write(iout,'(a,3(2x,f16.9))')' Lambda = ',(lambda(idir),idir = 1, 3)
522      write(iout,'(a,e16.9)')' Value of the energy functional E_1 = ',e1
523      write(iout,*)
524      write(iout,*)'Difference between actual value of the Polarization (C/m^2)'
525      write(iout,*)'and the target value:'
526    end if
527    diffpol(:) = (ptot_cart(:) - targetpol(:))*e_Cb/((Bohr_Ang*1.0d-10)**2)
528    if (iwrite) write(iout,'(3(3x,f16.9))')(diffpol(idir),idir = 1, 3)
529 
530    if (relaxat == 1) then
531 !    Compute the forces induced on the atoms by the electric field
532 !    The strength of the field is determined by lambda
533      felfd(:,:) = zero
534      do iatom = 1, natom
535        do idir = 1, 3
536          do jdir = 1, 3
537            index = idir +3*((iatom - 1) + mpert*((jdir - 1) + 3*(natom + 1)))
538            felfd(idir,iatom) = felfd(idir,iatom) - lambda(jdir)*blkval(1,index)/ucvol
539          end do
540        end do
541      end do
542 
543 !    Compute remaining forces and write them out
544 
545      fdiff(:,:) = fcart(:,:) - felfd(:,:)
546      if (iwrite) then
547        write(iout,*)
548        write(iout,*)'Difference between the Hellmann-Feynman forces'
549        write(iout,*)'and the forces induced by the electric field'
550        write(iout,*)'(cartesian coordinates, hartree/bohr)'
551      end if
552      fmax = zero
553      do iatom = 1, natom
554        if (iwrite) write(iout,'(3(3x,es16.9))')(fdiff(idir,iatom),idir = 1, 3)
555        do idir = 1, 3
556          if (abs(fdiff(idir,iatom)) > fmax) fmax = abs(fdiff(idir,iatom))
557        end do
558      end do
559 
560      if (iwrite) then
561        write(iout,'(a,3x,es16.9)')' fmax = ',fmax
562        write(iout,*)
563        write(iout,*)'Change of cartesian coordinates (delta_xcart):'
564        do iatom = 1, natom
565          write(iout,'(5x,i3,3(2x,f16.9))')iatom,(delta_xcart(idir,iatom),idir = 1, 3)
566        end do
567        write(iout,*)
568        write(iout,*)'New cartesian coordinates (xcart_new):'
569        write(iout,*)'  xcart'
570        do iatom = 1, natom
571          write(iout,'(3(3x,d22.14))')(xcart_new(idir,iatom),idir = 1, 3)
572        end do
573        write(iout,*)
574        write(iout,*)'New reduced coordinates (xred_new):'
575        write(iout,*)'  xred'
576        do iatom = 1, natom
577          write(iout,'(3(3x,d22.14))')(xred_new(idir,iatom),idir = 1, 3)
578        end do
579      end if
580 
581    end if         ! relaxat == 1
582 
583    if (relaxstr == 1) then
584 
585 !    Compute the stresses induced by the electric field
586      sigelfd(:) = zero
587      istrain = 0
588      do ipert = 1, 2
589        jpert = natom + 2 + ipert
590        do idir = 1, 3
591          istrain = istrain + 1
592          do jdir = 1, 3
593            index = idir +3*((jpert - 1) + mpert*((jdir - 1) + 3*(natom + 1)))
594            sigelfd(istrain) = sigelfd(istrain) + lambda(jdir)*blkval(1,index)
595          end do
596          sigelfd(istrain) = sigelfd(istrain)/ucvol
597        end do
598      end do
599 
600 !    Compute the remaining stresses and write them out
601      diffsig(:) = strten(:) - sigelfd(:)
602      sigmax = zero
603      do istrain = 1, 6
604        if (abs(diffsig(istrain)) > sigmax) sigmax = abs(diffsig(istrain))
605      end do
606      if (iwrite) then
607        write(iout,*)
608        write(iout,*)'Difference between the Hellmann-Feynman stresses'
609        write(iout,*)'and the stresses induced by the electric field'
610        write(iout,*)'(cartesian coordinates, hartree/bohr^3)'
611        write(iout,'(2x,a,f16.9,5x,a,f16.9)')'diffsig(1) = ',diffsig(1),'diffsig(4) = ',diffsig(4)
612        write(iout,'(2x,a,f16.9,5x,a,f16.9)')'diffsig(2) = ',diffsig(2),'diffsig(5) = ',diffsig(5)
613        write(iout,'(2x,a,f16.9,5x,a,f16.9)')'diffsig(3) = ',diffsig(3),'diffsig(6) = ',diffsig(6)
614        write(iout,'(a,3x,es16.9)')' sigmax = ',sigmax
615        write(iout,*)
616        write(iout,*)'Induced strain (delta_eta):'
617        write(iout,'(2x,a,f16.9,5x,a,f16.9)')'delta_eta(1) = ',delta_eta(1),'delta_eta(4) = ',delta_eta(4)
618        write(iout,'(2x,a,f16.9,5x,a,f16.9)')'delta_eta(2) = ',delta_eta(2),'delta_eta(5) = ',delta_eta(5)
619        write(iout,'(2x,a,f16.9,5x,a,f16.9)')'delta_eta(3) = ',delta_eta(3),'delta_eta(6) = ',delta_eta(6)
620        write(iout,*)
621        write(iout,*)'New lattice constants (acell_new):'
622        write(iout,*)'  acell'
623        write(iout,'(3(2x,d22.14))')(acell_new(idir),idir = 1, 3)
624        write(iout,*)
625        write(iout,*)'New primitive vectors (rprim_new):'
626        write(iout,*)'  rprim'
627        write(iout,'(3(2x,d22.14))')(rprim(idir,1),idir = 1, 3)
628        write(iout,'(3(2x,d22.14))')(rprim(idir,2),idir = 1, 3)
629        write(iout,'(3(2x,d22.14))')(rprim(idir,3),idir = 1, 3)
630      end if
631    end if         ! relaxstr /= 0
632 
633  end if    !  (relaxat /= 0).or.(relaxstr /= 0)
634 
635  ABI_DEALLOCATE(cfac)
636  ABI_DEALLOCATE(fdiff)
637  ABI_DEALLOCATE(felfd)
638  ABI_DEALLOCATE(delta)
639  ABI_DEALLOCATE(fcart)
640  ABI_DEALLOCATE(fcmat)
641  ABI_DEALLOCATE(ifcmat)
642  ABI_DEALLOCATE(ipvt)
643  ABI_DEALLOCATE(rfpert)
644  ABI_DEALLOCATE(vec)
645  ABI_DEALLOCATE(zgwork)
646  ABI_DEALLOCATE(irelaxat)
647 
648 end subroutine relaxpol