TABLE OF CONTENTS


ABINIT/m_relaxpol [ Modules ]

[ Top ] [ Modules ]

NAME

  m_relaxpol

FUNCTION

COPYRIGHT

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

SOURCE

15 #if defined HAVE_CONFIG_H
16 #include "config.h"
17 #endif
18 
19 #include "abi_common.h"
20 
21 module m_relaxpol
22 
23  use defs_basis
24  use m_abicore
25  use m_errors
26 
27  use m_fstrings,  only : sjoin, itoa
28  use m_symtk,     only : matr3inv
29  use m_berrytk,   only : polcart
30  use m_hide_lapack,   only : dzgedi, dzgefa
31  use m_geometry,  only : xcart2xred
32  use m_dynmat,    only : symdyma
33  use m_crystal,   only : crystal_t
34 
35  implicit none
36 
37  private

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) [[cite:Sai2002]].

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
 gred(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.

SOURCE

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