TABLE OF CONTENTS


ABINIT/m_relaxpol [ Modules ]

[ Top ] [ Modules ]

NAME

  m_relaxpol

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_relaxpol
27 
28  use defs_basis
29  use m_abicore
30  use m_errors
31 
32  use m_fstrings,  only : sjoin, itoa
33  use m_symtk,     only : matr3inv
34  use m_berrytk,   only : polcart
35  use m_hide_lapack,   only : dzgedi, dzgefa
36  use m_geometry,  only : xcart2xred
37  use m_dynmat,    only : symdyma
38  use m_crystal,   only : crystal_t
39 
40 
41  implicit none
42 
43  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
 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

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