TABLE OF CONTENTS
ABINIT/dfpt_gatherdy [ Functions ]
NAME
dfpt_gatherdy
FUNCTION
Sum (gather) the different parts of the 2nd-order matrix, to get the matrix of second-order derivatives (d2matr) Then, generates the dynamical matrix, not including the masses, but the correct non-cartesian coordinates ( => d2cart)
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (XG, DRH, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
becfrnl(3,natom,3*pawbec)=NL frozen contribution to Born Effective Charges (PAW only) berryopt=option for berry phase treatment blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical matrix has been calculated ; 0 otherwise ) dyew(2,3,natom,3,natom)=Ewald part of the dyn.matrix dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag)=frozen wf part of the dyn.matrix (except xc1) dyfrx1(2,3,natom,3,natom)=xc core correction (1) part of the frozen-wf part of the dynamical matrix. dyfr_cplex=1 if dyfrnl is real, 2 if it is complex dyfr_nondiag=1 if dyfrwf is non diagonal with respect to atoms; 0 otherwise dyvdw(2,3,natom,3,natom*usevdw)=vdw DFT-D part of the dynamical matrix d2bbb(2,3,3,mpert,mband,mband*prtbbb)=band by band decomposition of some second order derivatives d2nfr(2,3,mpert,3,mpert)=non-frozen wf part of the 2nd-order matr eltcore(6,6)=core contribution to the elastic tensor elteew(6+3*natom,6)=Ewald contribution to the elastic tsenor eltfrhar(6,6)=hartree contribution to the elastic tensor eltfrkin(6,6)=kinetic contribution to the elastic tensor eltfrloc(6+3*natom,6)=local psp contribution to the elastic tensor eltfrnl(6+3*natom,6)=non-local psp contribution to the elastic tensor eltfrxc(6+3*natom,6)=exchange-correlation contribution to the elastic tensor eltvdw(6+3*natom,6*usevdw)=vdw DFT-D part of the elastic tensor gprimd(3,3)=basis vector in the reciprocal space mband=maximum number of bands mpert =maximum number of ipert natom=number of atoms in unit cell ntypat=number of atom types outd2=option for the output of the 2nd-order matrix : if outd2=1, non-stationary part if outd2=2, stationary part. pawbec= flag for the computation of frozen part of Born Effective Charges (PAW only) prtbbb=if 1, print the band-by-band decomposition, otherwise, prtbbb=0 rfasr= (0=> no acoustic sum rule [asr] imposed), (1=> asr is imposed, in the democratic way for the effective charges), (2=> asr is imposed, in the aristocratic way for the effective charges) rfpert(mpert)=define the perturbations rprimd(3,3)=dimensional primitive translations (bohr) typat(natom)=integer label of each type of atom (1,2,...) ucvol=unit cell volume usevdw= flag set to 1 if vdw DFT-D semi-empirical potential is in use zion(ntypat)=charge corresponding to the atom type
OUTPUT
carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian 2DTE matrix has been calculated correctly ; 0 otherwise ) d2cart(2,3,mpert,3,mpert)= dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb)= band by band decomposition of Born effective charges (calculated from phonon-type perturbation) in cartesian coordinates d2matr(2,3,mpert,3,mpert)=2nd-order matrix (masses non included, no cartesian coordinates : simply second derivatives)
PARENTS
respfn
CHILDREN
asria_calc,asria_corr,cart29,cart39,chneu9
SOURCE
84 #if defined HAVE_CONFIG_H 85 #include "config.h" 86 #endif 87 88 #include "abi_common.h" 89 90 91 subroutine dfpt_gatherdy(becfrnl,berryopt,blkflg,carflg,dyew,dyfrwf,dyfrx1,& 92 & dyfr_cplex,dyfr_nondiag,dyvdw,d2bbb,d2cart,d2cart_bbb,d2matr,d2nfr,& 93 & eltcore,elteew,eltfrhar,eltfrkin,eltfrloc,eltfrnl,eltfrxc,eltvdw,& 94 & gprimd,mband,mpert,natom,ntypat,outd2,pawbec,pawpiezo,piezofrnl,prtbbb,& 95 & rfasr,rfpert,rprimd,typat,ucvol,usevdw,zion) 96 97 use defs_basis 98 use m_errors 99 use m_profiling_abi 100 101 use m_dynmat, only : asria_calc,asria_corr,cart29,cart39,chneu9 102 103 !This section has been created automatically by the script Abilint (TD). 104 !Do not modify the following lines by hand. 105 #undef ABI_FUNC 106 #define ABI_FUNC 'dfpt_gatherdy' 107 !End of the abilint section 108 109 implicit none 110 111 !Arguments ------------------------------- 112 !scalars 113 integer,intent(in) :: berryopt,dyfr_cplex,dyfr_nondiag,mband,mpert,natom,ntypat,outd2 114 integer,intent(in) :: pawbec,pawpiezo,prtbbb,rfasr,usevdw 115 real(dp),intent(in) :: ucvol 116 !arrays 117 integer,intent(in) :: rfpert(mpert),typat(natom) 118 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 119 integer,intent(out) :: carflg(3,mpert,3,mpert) 120 real(dp),intent(in) :: becfrnl(3,natom,3*pawbec) 121 real(dp),intent(in) :: d2bbb(2,3,3,mpert,mband,mband*prtbbb) 122 real(dp),intent(in) :: d2nfr(2,3,mpert,3,mpert),dyew(2,3,natom,3,natom) 123 real(dp),intent(in) :: dyfrwf(dyfr_cplex,3,3,natom,1+(natom-1)*dyfr_nondiag) 124 real(dp),intent(in) :: dyfrx1(2,3,natom,3,natom),dyvdw(2,3,natom,3,natom*usevdw) 125 real(dp),intent(in) :: eltcore(6,6),elteew(6+3*natom,6) 126 real(dp),intent(in) :: eltfrhar(6,6),eltfrkin(6,6),eltfrloc(6+3*natom,6) 127 real(dp),intent(in) :: eltfrnl(6+3*natom,6),eltfrxc(6+3*natom,6) 128 real(dp),intent(in) :: eltvdw(6+3*natom,6*usevdw),gprimd(3,3) 129 real(dp),intent(in) :: piezofrnl(6,3*pawpiezo),rprimd(3,3),zion(ntypat) 130 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert) 131 real(dp),intent(out) :: d2cart_bbb(2,3,3,mpert,mband,mband*prtbbb) 132 real(dp),intent(out) :: d2matr(2,3,mpert,3,mpert) 133 134 !Local variables ------------------------- 135 !scalars 136 integer :: chneut,iband,iblok,idir,idir1,idir2,ii,ipert,ipert1,ipert2 137 integer :: jj,nblok,selectz 138 character(len=500) :: message 139 !arrays 140 integer :: flg1(3),flg2(3) 141 real(dp) :: vec1(3),vec2(3) 142 ! real(dp) :: ter(3,3) ! this variable appears commented out below 143 real(dp),allocatable :: d2tmp(:,:,:,:,:),d2work(:,:,:,:,:),elfrtot(:,:) 144 145 ! ********************************************************************* 146 147 148 !DEBUG 149 !write(std_out,*)' dfpt_gatherdy : enter ' 150 !write(std_out,*)' outd2,mpert =',outd2,mpert 151 !write(std_out,*)' blkflg(:,natom+2,:,natom+2)=',blkflg(:,natom+2,:,natom+2) 152 !ENDDEBUG 153 154 if(outd2/=3)then 155 156 ! Initialise the 2nd-derivative matrix 157 d2matr(:,:,:,:,:)=0.0_dp 158 159 ! Add the non-frozen-part, the 160 ! Ewald part and the xc1 part of the frozen-wf part 161 ! Add the vdw part (if any) 162 do ipert2=1,mpert 163 do idir2=1,3 164 do ipert1=1,mpert 165 do idir1=1,3 166 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 167 do ii=1,2 168 d2matr(ii,idir1,ipert1,idir2,ipert2)=& 169 & d2nfr(ii,idir1,ipert1,idir2,ipert2) 170 if(ipert1<=natom .and. ipert2<=natom) then 171 d2matr(ii,idir1,ipert1,idir2,ipert2)=& 172 & d2matr(ii,idir1,ipert1,idir2,ipert2)+& 173 & dyew(ii,idir1,ipert1,idir2,ipert2) +& 174 & dyfrx1(ii,idir1,ipert1,idir2,ipert2) 175 if (usevdw==1) then 176 d2matr(ii,idir1,ipert1,idir2,ipert2)=& 177 & d2matr(ii,idir1,ipert1,idir2,ipert2)+& 178 & dyvdw(ii,idir1,ipert1,idir2,ipert2) 179 end if 180 end if 181 end do 182 end if 183 end do 184 end do 185 end do 186 end do 187 188 ! Add the frozen-wavefunction part 189 if (dyfr_nondiag==0) then 190 do ipert2=1,natom 191 do idir2=1,3 192 do idir1=1,3 193 if( blkflg(idir1,ipert2,idir2,ipert2)==1 ) then 194 d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)=& 195 & d2matr(1:dyfr_cplex,idir1,ipert2,idir2,ipert2)& 196 & +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert2,1) 197 end if 198 end do 199 end do 200 end do 201 else 202 do ipert2=1,natom 203 do ipert1=1,natom 204 do idir2=1,3 205 do idir1=1,3 206 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 207 d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)=& 208 & d2matr(1:dyfr_cplex,idir1,ipert1,idir2,ipert2)& 209 & +dyfrwf(1:dyfr_cplex,idir1,idir2,ipert1,ipert2) 210 end if 211 end do 212 end do 213 end do 214 end do 215 end if 216 217 ! Add the frozen-wavefunction part of Born Effective Charges 218 if (pawbec==1) then 219 ipert2=natom+2 220 do idir2=1,3 ! Direction of electric field 221 do ipert1=1,natom ! Atom 222 do idir1=1,3 ! Direction of atom 223 if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 224 d2matr(1,idir1,ipert1,idir2,ipert2)=& 225 & d2matr(1,idir1,ipert1,idir2,ipert2)+becfrnl(idir1,ipert1,idir2) 226 end if 227 if(blkflg(idir2,ipert2,idir1,ipert1)==1 ) then 228 d2matr(1,idir2,ipert2,idir1,ipert1)=& 229 & d2matr(1,idir2,ipert2,idir1,ipert1)+becfrnl(idir1,ipert1,idir2) 230 end if 231 end do 232 end do 233 end do 234 end if 235 236 ! Section for piezoelectric tensor (from electric field response only for PAW) 237 if(rfpert(natom+2)==1.and.pawpiezo==1) then 238 ipert2=natom+2 239 do idir2=1,3 ! Direction of electric field 240 do ipert1=natom+3,natom+4 ! Strain 241 do idir1=1,3 242 ii=idir1+3*(ipert1-natom-3) 243 if(blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 244 d2matr(1,idir1,ipert1,idir2,ipert2)=& 245 & d2nfr(1,idir1,ipert1,idir2,ipert2)+piezofrnl(ii,idir2) 246 end if 247 end do 248 end do 249 end do 250 end if 251 252 ! Section for strain perturbation 253 if(rfpert(natom+3)==1 .or. rfpert(natom+4)==1) then 254 ! Make sure relevant columns of output are nulled 255 d2matr(:,:,:,:,natom+3:natom+4)=0.0_dp 256 ! Accumulate all frozen parts of the elastic tensor 257 ABI_ALLOCATE(elfrtot,(6+3*natom,6)) 258 elfrtot(:,:)=elteew(:,:)+eltfrloc(:,:)+eltfrnl(:,:)+eltfrxc(:,:) 259 elfrtot(1:6,1:6)=elfrtot(1:6,1:6)+eltcore(:,:)+eltfrhar(:,:)+eltfrkin(:,:) 260 if (usevdw==1) elfrtot(:,:)=elfrtot(:,:)+eltvdw(:,:) 261 262 do ipert2=natom+3,natom+4 263 do idir2=1,3 264 ! Internal strain components first 265 do ipert1=1,natom 266 do idir1=1,3 267 if( blkflg(1,ipert1,idir2,ipert2)==1 ) then 268 ii=idir1+6+3*(ipert1-1) 269 jj=idir2+3*(ipert2-natom-3) 270 d2matr(1,idir1,ipert1,idir2,ipert2)=& 271 & d2nfr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj) 272 end if 273 end do 274 end do 275 ! Now, electric field - strain mixed derivative (piezoelectric tensor) 276 ipert1=natom+2 277 do idir1=1,3 278 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 279 d2matr(1,idir1,ipert1,idir2,ipert2)=& 280 & d2nfr(1,idir1,ipert1,idir2,ipert2) 281 if (pawpiezo==1) then 282 ii=idir2+3*(ipert2-natom-3) 283 d2matr(1,idir1,ipert1,idir2,ipert2)=& 284 & d2matr(1,idir1,ipert1,idir2,ipert2)+piezofrnl(ii,idir1) 285 end if 286 end if 287 end do 288 ! Now, strain-strain 2nd derivatives 289 do ipert1=natom+3,natom+4 290 do idir1=1,3 291 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 292 ii=idir1+3*(ipert1-natom-3) 293 jj=idir2+3*(ipert2-natom-3) 294 d2matr(1,idir1,ipert1,idir2,ipert2)=& 295 & d2nfr(1,idir1,ipert1,idir2,ipert2)+elfrtot(ii,jj) 296 end if 297 end do 298 end do 299 end do 300 end do 301 ABI_DEALLOCATE(elfrtot) 302 end if 303 ! End section for strain perturbation 304 305 ! The second-order matrix has been computed. 306 307 ! Filter now components smaller in absolute value than 1.0d-20, 308 ! for automatic testing reasons 309 do ipert2=1,mpert 310 do idir2=1,3 311 do ipert1=1,mpert 312 do idir1=1,3 313 if( blkflg(idir1,ipert1,idir2,ipert2)==1 ) then 314 do ii=1,2 315 if( d2matr(ii,idir1,ipert1,idir2,ipert2)**2 < 1.0d-40)then 316 d2matr(ii,idir1,ipert1,idir2,ipert2)=zero 317 end if 318 end do 319 end if 320 end do 321 end do 322 end do 323 end do 324 325 ! DEBUG 326 ! write(std_out,*)' d2matr ' 327 ! ipert2=natom+2 328 ! do idir2=1,3 329 ! ipert1=natom+2 330 ! do idir1=1,3 331 ! write(std_out,'(4i4,2es16.6)' )idir1,ipert1,idir2,ipert2,& 332 ! & d2matr(1,idir1,ipert1,idir2,ipert2),& 333 ! & d2matr(2,idir1,ipert1,idir2,ipert2) 334 ! end do 335 ! end do 336 ! ENDDEBUG 337 338 ! Cartesian coordinates transformation 339 iblok=1 ; nblok=1 340 341 ! In the case of finite electric field, the convention for the 342 ! direction of the electric field perturbation was NOT the usual convention ... 343 ! So, there is a transformation to the usual conventions 344 ! to be done first ... 345 if((berryopt==4 .or. berryopt==6 .or. berryopt==7 .or. berryopt==14 .or. berryopt==16 .or. berryopt==17 ) & 346 & .and. minval(abs(blkflg(:,natom+2,:,natom+2)))/=0)then !!HONG need to check for fixed D and E calculation 347 if(minval(abs(blkflg(:,natom+2,:,natom+2)-1))/=0)then 348 write(message,'(5a)')& 349 & ' In case of finite electric field, and electric field perturbation,',ch10,& 350 & ' the three directions for the perturbations must be treated.',ch10,& 351 & ' Action : set idir to 1 1 1, or forget about finite electric field.' 352 MSG_ERROR(message) 353 end if 354 do ipert=1,mpert 355 do idir=1,3 356 do ii=1,2 357 vec1(:)=d2matr(ii,idir,ipert,:,natom+2) 358 flg1(:)=blkflg(idir,ipert,:,natom+2) 359 call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2) 360 d2matr(ii,idir,ipert,:,natom+2)=vec2(:)*two_pi 361 blkflg(idir,ipert,:,natom+2)=flg2(:) 362 end do 363 end do 364 end do 365 do ipert=1,mpert 366 do idir=1,3 367 do ii=1,2 368 vec1(:)=d2matr(ii,:,natom+2,idir,ipert) 369 flg1(:)=blkflg(:,natom+2,idir,ipert) 370 call cart39(flg1,flg2,gprimd,1,1,rprimd,vec1,vec2) 371 d2matr(ii,:,natom+2,idir,ipert)=vec2(:)*two_pi 372 blkflg(:,natom+2,idir,ipert)=flg2(:) 373 end do 374 end do 375 end do 376 ! Also to be done, a change of sign, that I do not understand (XG071110) 377 ! Perhaps due to d/dk replacing id/dk ? ! 378 d2matr(:,:,natom+2,:,natom+2)=-d2matr(:,:,natom+2,:,natom+2) 379 end if 380 381 call cart29(blkflg,d2matr,carflg,d2cart,& 382 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion) 383 384 ! Band by band decomposition of the Born effective charges 385 if(prtbbb==1)then 386 ABI_ALLOCATE(d2work,(2,3,mpert,3,mpert)) 387 ABI_ALLOCATE(d2tmp,(2,3,mpert,3,mpert)) 388 do iband=1,mband 389 d2work(:,:,:,:,:)=0.0_dp 390 d2tmp(:,:,:,:,:)=0.0_dp 391 d2work(:,:,natom+2,:,:) = d2bbb(:,:,:,:,iband,iband) 392 call cart29(blkflg,d2work,carflg,d2tmp,& 393 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion) 394 395 ! Remove the ionic part 396 do ipert1=1,natom 397 do idir1=1,3 398 d2tmp(1,idir1,natom+2,idir1,ipert1) = & 399 & d2tmp(1,idir1,natom+2,idir1,ipert1) - zion(typat(ipert1)) 400 end do 401 end do 402 403 d2cart_bbb(:,:,:,:,iband,iband)=d2tmp(:,:,natom+2,:,:) 404 405 end do 406 ABI_DEALLOCATE(d2tmp) 407 ABI_DEALLOCATE(d2work) 408 end if ! prtbbb==1 409 410 ! 411 ! Now, the cartesian elements are ready for output 412 ! carflg give the information on their correctness 413 end if 414 415 ! Imposition of the ASR on the analytical part of the DynMat 416 ! Assume that if rfasr/=0, the whole cartesian matrix is correct 417 if(rfasr/=0)then 418 419 ABI_ALLOCATE(d2work,(2,3,mpert,3,mpert)) 420 call asria_calc(rfasr,d2work,d2cart,mpert,natom) 421 ! The following line imposes ASR: 422 call asria_corr(rfasr,d2work,d2cart,mpert,natom) 423 424 ABI_DEALLOCATE(d2work) 425 426 ! Imposition of the ASR on the effective charges. 427 if(rfpert(natom+2)==1)then 428 chneut=rfasr 429 selectz=0 430 call chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion) 431 end if 432 433 end if 434 435 !Additional operations on cartesian strain derivatives 436 if(rfpert(natom+3)==1 .or. rfpert(natom+4)==1) then 437 ! Impose zero-net-force condition on internal strain tensor 438 do ipert2=natom+3,natom+4 439 do idir2=1,3 440 vec1(:)=0.0_dp 441 do ipert1=1,natom 442 do idir1=1,3 443 if(carflg(idir1,ipert1,idir2,ipert2)==1) then 444 vec1(idir1)=vec1(idir1)+d2cart(1,idir1,ipert1,idir2,ipert2) 445 end if 446 end do 447 end do 448 vec1(:)=vec1(:)/dble(natom) 449 do ipert1=1,natom 450 do idir1=1,3 451 if(carflg(idir1,ipert1,idir2,ipert2)==1) then 452 ! Note minus sign to convert gradients to forces 453 d2cart(1,idir1,ipert1,idir2,ipert2)=& 454 & -(d2cart(1,idir1,ipert1,idir2,ipert2)-vec1(idir1)) 455 end if 456 end do 457 end do 458 end do 459 end do 460 ! Divide strain 2nd deriviative by ucvol to give elastic tensor 461 do ipert2=natom+3,natom+4 462 do idir2=1,3 463 do ipert1=natom+3,natom+4 464 do idir1=1,3 465 if(carflg(idir1,ipert1,idir2,ipert2)==1) then 466 d2cart(1,idir1,ipert1,idir2,ipert2)=& 467 & d2cart(1,idir1,ipert1,idir2,ipert2)/ucvol 468 end if 469 end do 470 end do 471 end do 472 end do 473 end if 474 475 !calculate Born effective charges from electric field perturbation 476 !do ipert1=1,natom 477 !ter(:,:)=zero 478 !do idir1=1,3 479 !do ii=1,3 480 !do jj=1,3 481 !if(abs(gprimd(idir1,ii))>1.0d-10)then 482 !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,idir1,natom+2,jj,ipert1)*gprimd(jj,ii) 483 !endif 484 !enddo 485 !enddo 486 !add zion to bec 487 !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1)) 488 !enddo 489 !d2cart(1,:,ipert1,:,natom+2)=ter(:,:) 490 !enddo 491 !carflg(:,1:natom,:,natom+2)=1 492 493 !Born effective charges from phonon perturbation 494 !do ipert1=1,natom 495 !ter(:,:)=zero 496 !do idir1=1,3 497 !do ii=1,3 498 !do jj=1,3 499 !if(abs(gprimd(idir1,ii))>1.0d-10)then 500 !ter(idir1,ii)=ter(idir1,ii)+ d2nfr(1,jj,ipert1,idir1,natom+2)*gprimd(jj,ii) 501 !endif 502 !enddo 503 !enddo 504 !! add zion to bec 505 !ter(idir1,idir1)=ter(idir1,idir1)+zion(typat(ipert1)) 506 !enddo 507 !d2cart(1,:,natom+2,:,ipert1)=ter(:,:) 508 !enddo 509 !carflg(:,natom+2,:,1:natom)=1 510 511 512 !!Dielectric constant 513 !do ii=1,3 514 !do jj=1,3 515 !ter(ii,jj)=d2nfr(1,ii,natom+2,jj,natom+2) 516 !end do 517 !end do 518 !ter(:,:)=pi*four*ter(:,:)/ucvol 519 ! 520 !do ii=1,3 521 !ter(ii,ii)=ter(ii,ii)+one 522 !end do 523 !d2cart(1,:,natom+2,:,natom+2)=ter(:,:) 524 !carflg(:,natom+2,:,1:natom+2)=1 525 526 !DEBUG 527 !Debugging, but needed to make the stuff work on the IBM Dirac ? ! 528 !write(std_out,*)' d2cart ' 529 !ipert2=natom+2 530 !do idir2=1,3 531 !ipert1=natom+2 532 !do idir1=1,3 533 !write(std_out,'(5i4,2d20.10)' )idir1,ipert1,idir2,ipert2,& 534 !& carflg(idir1,ipert1,idir2,ipert2),& 535 !& d2cart(1,idir1,ipert1,idir2,ipert2),& 536 !& d2cart(2,idir1,ipert1,idir2,ipert2) 537 !end do 538 !end do 539 !ENDDEBUG 540 541 !DEBUG 542 !write(std_out,*)' dfpt_gatherdy : exit ' 543 !ENDDEBUG 544 545 end subroutine dfpt_gatherdy