TABLE OF CONTENTS
ABINIT/m_relaxpol [ 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 ]
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