TABLE OF CONTENTS
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).
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