TABLE OF CONTENTS
ABINIT/nonlinear [ Functions ]
NAME
nonlinear
FUNCTION
Primary routine for conducting DFT calculations of non linear response functions.
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (MVeithen, MB) 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
codvsn = code version dtfil <type(datafiles_type)> = variables related to files dtset <type(dataset_type)> = all input variables for this dataset etotal = new total energy (no meaning at output) iexit= exit flag mband = maximum number of bands mgfft = maximum single fft dimension mkmem = maximum number of k points which can fit in core memory mpi_enreg=informations about MPI pnarallelization mpw = maximum number of planewaves in basis sphere (large number) natom = number of atoms in unit cell nfft = (effective) number of FFT grid points (for this processor) nkpt = number of k points nspden = number of spin-density components nspinor = number of spinorial components of the wavefunctions nsppol = number of channels for spin-polarization (1 or 2) nsym = number of symmetry elements in space group occ(mband*nkpt*nsppol) = occupation number for each band and k xred(3,natom) = reduced atomic coordinates
OUTPUT
npwtot(nkpt) = total number of plane waves at each k point
SIDE EFFECTS
pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)> = variables related to pseudopotentials
PARENTS
driver
CHILDREN
d3sym,ddb_hdr_free,ddb_hdr_init,ddb_hdr_open_write,dfptnl_doutput dfptnl_loop,ebands_free,fourdp,getcut,getkgrid,getshell,hartre,hdr_free hdr_init,hdr_update,initmv,inwffil,kpgio,mkcore,nlopt,pspini,read_rhor rhotoxc,setsym,setup1,status,symmetrize_xred,sytens,timab,wffclose wrtout,xcdata_init
SOURCE
59 #if defined HAVE_CONFIG_H 60 #include "config.h" 61 #endif 62 63 #include "abi_common.h" 64 65 subroutine nonlinear(codvsn,dtfil,dtset,etotal,iexit,& 66 & mband,mgfft,mkmem,mpi_enreg,mpw,natom,nfft,nkpt,npwtot,nspden,& 67 & nspinor,nsppol,nsym,occ,pawrad,pawtab,psps,xred) 68 69 use defs_basis 70 use defs_datatypes 71 use defs_abitypes 72 use defs_wvltypes 73 use m_wffile 74 use m_errors 75 use m_profiling_abi 76 use m_xmpi 77 use m_hdr 78 use m_ebands 79 use m_xcdata 80 81 use m_dynmat, only : d3sym, sytens 82 use m_ddb, only : nlopt, DDB_VERSION 83 use m_ddb_hdr, only : ddb_hdr_type, ddb_hdr_init, ddb_hdr_free, ddb_hdr_open_write 84 use m_ioarr, only : read_rhor 85 use m_pawrad, only : pawrad_type 86 use m_pawtab, only : pawtab_type 87 use m_pawrhoij, only : pawrhoij_type 88 use m_kg, only : getcut, kpgio 89 90 !This section has been created automatically by the script Abilint (TD). 91 !Do not modify the following lines by hand. 92 #undef ABI_FUNC 93 #define ABI_FUNC 'nonlinear' 94 use interfaces_14_hidewrite 95 use interfaces_18_timing 96 use interfaces_32_util 97 use interfaces_41_geometry 98 use interfaces_53_ffts 99 use interfaces_56_recipspace 100 use interfaces_56_xc 101 use interfaces_64_psp 102 use interfaces_67_common 103 use interfaces_72_response 104 use interfaces_79_seqpar_mpi 105 use interfaces_95_drive, except_this_one => nonlinear 106 !End of the abilint section 107 108 implicit none 109 110 !Arguments ------------------------------------ 111 !scalars 112 integer,intent(in) :: iexit,mband,mgfft,mkmem,mpw,nfft 113 integer,intent(in) :: natom,nkpt,nspden,nspinor,nsppol,nsym 114 logical :: non_magnetic_xc 115 real(dp),intent(inout) :: etotal 116 character(len=6),intent(in) :: codvsn 117 type(MPI_type),intent(inout) :: mpi_enreg 118 type(datafiles_type),intent(in) :: dtfil 119 type(dataset_type),intent(inout) :: dtset 120 type(pseudopotential_type),intent(inout) :: psps 121 !arrays 122 integer,intent(out) :: npwtot(nkpt) 123 real(dp),intent(inout) :: occ(mband*nkpt*nsppol),xred(3,natom) 124 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat,psps%usepaw) 125 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat,psps%usepaw) 126 127 !Local variables------------------------------- 128 !scalars 129 integer,parameter :: level=50,response=1,cplex1=1 130 integer :: ask_accurate,bantot,dum_nshiftk,flag 131 integer :: formeig,gencond,gscase,i1dir,i1pert,i2dir,i2pert,i3dir 132 integer :: i3pert,ierr,ireadwf,mcg,mkmem_max,mpert,n1,n2,n3,n3xccc 133 integer :: nkpt3,nkxc,nk3xc,nneigh,option,optorth,rdwrpaw,comm_cell 134 real(dp) :: boxcut,ecore,ecut_eff,enxc,fermie,gsqcut,gsqcut_eff,gsqcut_eff2,gsqcutdg_eff 135 real(dp) :: rdum,residm,ucvol,vxcavg 136 character(len=500) :: message 137 character(len=fnlen) :: dscrpt 138 type(ebands_t) :: bstruct 139 type(hdr_type) :: hdr,hdr_den 140 type(ddb_hdr_type) :: ddb_hdr 141 type(wffile_type) :: wffgs,wfftgs 142 type(wvl_data) :: wvl 143 type(xcdata_type) :: xcdata 144 !arrays 145 integer :: dum_kptrlatt(3,3),dum_vacuum(3),perm(6) 146 integer,allocatable :: blkflg(:,:,:,:,:,:),carflg(:,:,:,:,:,:),cgindex(:,:) 147 integer,allocatable :: d3e_pert1(:),d3e_pert2(:),d3e_pert3(:) 148 integer,allocatable :: indsym(:,:,:),irrzon(:,:,:),kg(:,:),kneigh(:,:),kg_neigh(:,:,:) 149 integer,allocatable :: kptindex(:,:),npwarr(:),pwind(:,:,:),rfpert(:,:,:,:,:,:),symrec(:,:,:) 150 real(dp) :: dum_shiftk(3,210),dummy2(6),gmet(3,3),gprimd(3,3),k0(3) 151 real(dp) :: rmet(3,3),rprimd(3,3),strsxc(6),tsec(2) 152 real(dp),allocatable :: amass(:),cg(:,:),d3cart(:,:,:,:,:,:,:) 153 real(dp),allocatable :: d3lo(:,:,:,:,:,:,:),dum_kptns(:,:) 154 real(dp),allocatable :: dum_wtk(:),dyfrx2(:,:,:),eigen(:),grxc(:,:),k3xc(:,:) 155 real(dp),allocatable :: kpt3(:,:),kxc(:,:),mvwtk(:,:),phnons(:,:,:),rhog(:,:) 156 real(dp),allocatable :: rhor(:,:),vhartr(:),vxc(:,:),work(:),xccc3d(:) 157 type(pawrhoij_type),allocatable :: pawrhoij(:) 158 159 ! *********************************************************************** 160 161 call timab(501,1,tsec) 162 call status(0,dtfil%filstat,iexit,level,'enter ') 163 164 comm_cell = mpi_enreg%comm_cell 165 166 ! Initialise non_magnetic_xc for rhohxc 167 non_magnetic_xc=(dtset%usepawu==4).or.(dtset%usepawu==14) 168 169 !Check if the perturbations asked in the input file can be computed 170 171 if (((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert2_phon == 1)).or. & 172 & ((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert3_phon == 1)).or. & 173 & ((dtset%d3e_pert2_phon == 1).and.(dtset%d3e_pert3_phon == 1))) then 174 write(message,'(7a)')& 175 & 'You have asked for a third-order derivative with respect to',ch10,& 176 & '2 or more atomic displacements.',ch10,& 177 & 'This is not allowed yet.',ch10,& 178 & 'Action : change d3e_pert1_phon, d3e_pert2_phon or d3e_pert3_phon in your input file.' 179 MSG_ERROR(message) 180 end if 181 182 !Define the set of admitted perturbations taking into account 183 !the possible permutations 184 mpert=natom+6 185 ABI_ALLOCATE(blkflg,(3,mpert,3,mpert,3,mpert)) 186 ABI_ALLOCATE(carflg,(3,mpert,3,mpert,3,mpert)) 187 ABI_ALLOCATE(rfpert,(3,mpert,3,mpert,3,mpert)) 188 ABI_ALLOCATE(d3e_pert1,(mpert)) 189 ABI_ALLOCATE(d3e_pert2,(mpert)) 190 ABI_ALLOCATE(d3e_pert3,(mpert)) 191 ABI_ALLOCATE(d3lo,(2,3,mpert,3,mpert,3,mpert)) 192 ABI_ALLOCATE(d3cart,(2,3,mpert,3,mpert,3,mpert)) 193 blkflg(:,:,:,:,:,:) = 0 194 d3lo(:,:,:,:,:,:,:) = 0_dp 195 rfpert(:,:,:,:,:,:) = 0 196 d3e_pert1(:) = 0 ; d3e_pert2(:) = 0 ; d3e_pert3(:) = 0 197 198 if (dtset%d3e_pert1_phon==1) d3e_pert1(dtset%d3e_pert1_atpol(1):dtset%d3e_pert1_atpol(2))=1 199 if (dtset%d3e_pert2_phon==1) d3e_pert2(dtset%d3e_pert2_atpol(1):dtset%d3e_pert2_atpol(2))=1 200 if (dtset%d3e_pert3_phon==1) d3e_pert3(dtset%d3e_pert3_atpol(1):dtset%d3e_pert3_atpol(2))=1 201 if (dtset%d3e_pert1_elfd/=0) d3e_pert1(natom+2)=1 202 if (dtset%d3e_pert2_elfd/=0) d3e_pert2(natom+2)=1 203 if (dtset%d3e_pert3_elfd/=0) d3e_pert3(natom+2)=1 204 205 do i1pert = 1, mpert 206 do i1dir = 1, 3 207 do i2pert = 1, mpert 208 do i2dir = 1, 3 209 do i3pert = 1, mpert 210 do i3dir = 1, 3 211 perm(1) = & 212 & d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) & 213 & *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) & 214 & *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir) 215 perm(2) = & 216 & d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) & 217 & *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) & 218 & *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir) 219 perm(3) = & 220 & d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) & 221 & *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) & 222 & *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir) 223 perm(4) = & 224 & d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) & 225 & *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) & 226 & *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir) 227 perm(5) = & 228 & d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) & 229 & *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) & 230 & *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir) 231 perm(6) = & 232 & d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) & 233 & *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) & 234 & *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir) 235 if (sum(perm(:)) > 0) rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 236 end do 237 end do 238 end do 239 end do 240 end do 241 end do 242 243 !Determine the symmetrical perturbations 244 ABI_ALLOCATE(irrzon,(dtset%nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))) 245 ABI_ALLOCATE(phnons,(2,dtset%nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))) 246 ABI_ALLOCATE(indsym,(4,nsym,natom)) 247 ABI_ALLOCATE(symrec,(3,3,nsym)) 248 call status(0,dtfil%filstat,iexit,level,'call setsym ') 249 call setsym(indsym,irrzon,dtset%iscf,natom,& 250 & nfft,dtset%ngfft,nspden,nsppol,nsym,& 251 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred) 252 253 call status(0,dtfil%filstat,iexit,level,'call sym_xred ') 254 call symmetrize_xred(indsym,natom,nsym,dtset%symrel,dtset%tnons,xred) 255 call status(0,dtfil%filstat,iexit,level,'call sytens ') 256 call sytens(indsym,mpert,natom,nsym,rfpert,symrec,dtset%symrel) 257 258 write(message, '(a,a,a,a,a)' ) ch10, & 259 & ' The list of irreducible elements of the Raman and non-linear',& 260 & ch10,' optical susceptibility tensors is:',ch10 261 call wrtout(ab_out,message,'COLL') 262 call wrtout(std_out,message,'COLL') 263 264 write(message,'(12x,a)')& 265 & 'i1pert i1dir i2pert i2dir i3pert i3dir' 266 call wrtout(ab_out,message,'COLL') 267 call wrtout(std_out,message,'COLL') 268 n1 = 0 269 do i1pert = 1, natom + 2 270 do i1dir = 1, 3 271 do i2pert = 1, natom + 2 272 do i2dir = 1, 3 273 do i3pert = 1, natom + 2 274 do i3dir = 1,3 275 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 276 n1 = n1 + 1 277 write(message,'(2x,i4,a,6(5x,i3))') n1,')', & 278 & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir 279 call wrtout(ab_out,message,'COLL') 280 call wrtout(std_out,message,'COLL') 281 else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) then 282 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 283 end if 284 end do 285 end do 286 end do 287 end do 288 end do 289 end do 290 write(message,'(a,a)') ch10,ch10 291 call wrtout(ab_out,message,'COLL') 292 call wrtout(std_out,message,'COLL') 293 294 !if (dtset%paral_rf == -1) then 295 write(std_out,'(a)')"--- !IrredPerts" 296 write(std_out,'(a)')'# List of irreducible perturbations for nonlinear' 297 write(std_out,'(a)')'irred_perts:' 298 299 n1 = 0 300 do i1pert = 1, natom + 2 301 do i1dir = 1, 3 302 do i2pert = 1, natom + 2 303 do i2dir = 1, 3 304 do i3pert = 1, natom + 2 305 do i3dir = 1,3 306 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 307 n1 = n1 + 1 308 write(std_out,'(a,i0)')" - i1pert: ",i1pert 309 write(std_out,'(a,i0)')" i1dir: ",i1dir 310 write(std_out,'(a,i0)')" i2pert: ",i2pert 311 write(std_out,'(a,i0)')" i2dir: ",i2dir 312 write(std_out,'(a,i0)')" i3pert: ",i3pert 313 write(std_out,'(a,i0)')" i3dir: ",i3dir 314 end if 315 end do 316 end do 317 end do 318 end do 319 end do 320 end do 321 write(std_out,'(a)')"..." 322 !MSG_ERROR_NODUMP("aborting now") 323 !end if 324 325 !Set up for iterations 326 ecut_eff= (dtset%ecut) * (dtset%dilatmx) **2 327 ABI_ALLOCATE(amass,(natom)) 328 call status(0,dtfil%filstat,iexit,level,'call setup1 ') 329 call setup1(dtset%acell_orig(1:3,1),amass,dtset%amu_orig(:,1),bantot,dtset,& 330 & ecut_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcut_eff2,& 331 & natom,dtset%ngfft,dtset%ngfft,nkpt,nsppol,& 332 & response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw) 333 334 !Set up the basis sphere of planewaves 335 ABI_ALLOCATE(kg,(3,mpw*dtset%mk1mem)) 336 ABI_ALLOCATE(npwarr,(nkpt)) 337 call status(0,dtfil%filstat,iexit,level,'call kpgio ') 338 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg,& 339 & dtset%kptns,mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,& 340 & mpw,npwarr,npwtot,nsppol) 341 342 !Recompute first large sphere cut-off gsqcut, 343 !without taking into account dilatmx 344 k0(:)=0.0_dp 345 call status(0,dtfil%filstat,iexit,level,'call getcut ') 346 call getcut(boxcut,dtset%ecut,gmet,gsqcut,dtset%iboxcut,std_out,k0,dtset%ngfft) 347 348 !Open and read pseudopotential files 349 ecore = 0_dp 350 call status(0,dtfil%filstat,iexit,level,'call pspini ') 351 call pspini(dtset,dtfil,ecore,gencond,gsqcut_eff,gsqcutdg_eff,& 352 & pawrad,pawtab,psps,rprimd,comm_mpi=mpi_enreg%comm_cell) 353 354 !Initialize band structure datatype 355 bstruct = ebands_from_dtset(dtset, npwarr) 356 357 !Initialize header 358 gscase=0 359 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr) 360 361 !Update header, with evolving variables, when available 362 !Here, rprimd, xred and occ are available 363 residm=hdr%residm ; fermie=hdr%fermie 364 call hdr_update(hdr,bantot,etotal,fermie,residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1)) 365 366 !Clean band structure datatype (should use it more in the future !) 367 call ebands_free(bstruct) 368 369 !Read ground-state wavefunctions 370 mcg=dtset%mpw*dtset%nspinor*mband*dtset%mkmem*dtset%nsppol 371 ABI_ALLOCATE(cg,(2,mcg)) 372 ABI_ALLOCATE(eigen,(mband*dtset%nkpt*dtset%nsppol)) 373 optorth=1;if (psps%usepaw==1) optorth=0 374 ireadwf=1 ; formeig=0 375 eigen(:)=0_dp ; ask_accurate=1 376 call status(0,dtfil%filstat,iexit,level,'call inwffil ') 377 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen,dtset%exchn2n3d,& 378 & formeig,hdr,ireadwf,dtset%istwfk,kg,dtset%kptns,& 379 & dtset%localrdwf,mband,mcg,dtset%mkmem,mpi_enreg,mpw,& 380 & dtset%nband,dtset%ngfft,dtset%nkpt,& 381 & npwarr,dtset%nsppol,dtset%nsym,& 382 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 383 & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,& 384 & dtfil%fnamewffk,wvl) 385 386 if (ireadwf==1) then 387 call WffClose(wffgs,ierr) 388 end if 389 390 ABI_DEALLOCATE(eigen) 391 392 ABI_ALLOCATE(rhog,(2,nfft)) 393 ABI_ALLOCATE(rhor,(nfft,nspden)) 394 !Get the ground state charge density 395 396 if (dtset%getden /= 0 .or. dtset%irdden /= 0) then 397 rdwrpaw = 0 398 call status(0,dtfil%filstat,iexit,level,'call ioarr ') 399 400 call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfft, dtset%ngfft, rdwrpaw, & 401 mpi_enreg, rhor, hdr_den, pawrhoij, comm_cell, check_hdr=hdr) 402 call hdr_free(hdr_den) 403 404 ! Compute up+down rho(G) by fft 405 ABI_ALLOCATE(work,(nfft)) 406 work(:)=rhor(:,1) 407 call status(0,dtfil%filstat,iexit,level,'call fourdp ') 408 call fourdp(1,rhog,work,-1,mpi_enreg,nfft,dtset%ngfft,dtset%paral_kgb,0) 409 ABI_DEALLOCATE(work) 410 end if 411 412 !Compute core electron density xccc3d 413 n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3) 414 ABI_ALLOCATE(grxc,(3,natom)) 415 ABI_ALLOCATE(vxc,(nfft,nspden)) 416 ABI_ALLOCATE(vhartr,(nfft)) 417 n3xccc=0 418 if (psps%n1xccc/=0) n3xccc=nfft 419 ABI_ALLOCATE(xccc3d,(n3xccc)) 420 if (psps%n1xccc/=0) then 421 option=1 422 ABI_ALLOCATE(dyfrx2,(3,3,natom)) 423 call status(0,dtfil%filstat,iexit,level,'call mkcore ') 424 call mkcore(dummy2,dyfrx2,grxc,mpi_enreg,natom,nfft,nspden,psps%ntypat,& 425 & n1,psps%n1xccc,n2,n3,option,rprimd,dtset%typat,ucvol,vxc,& 426 & psps%xcccrc,psps%xccc1d,xccc3d,xred) 427 ABI_DEALLOCATE(dyfrx2) 428 end if 429 430 call hartre(1,gsqcut,psps%usepaw,mpi_enreg,nfft,dtset%ngfft,dtset%paral_kgb,rhog,rprimd,vhartr) 431 432 !Compute kxc (second- and third-order exchange-correlation kernel) 433 option=3 434 nkxc=2*nspden-1 ! LDA 435 if(dtset%xclevel==2.and.nspden==1) nkxc=7 ! non-polarized GGA 436 if(dtset%xclevel==2.and.nspden==2) nkxc=19 ! polarized GGA 437 nk3xc=3*nspden-2 438 ABI_ALLOCATE(kxc,(nfft,nkxc)) 439 ABI_ALLOCATE(k3xc,(nfft,nk3xc)) 440 441 call status(0,dtfil%filstat,iexit,level,'call rhotoxc ') 442 ABI_ALLOCATE(work,(0)) 443 call xcdata_init(xcdata,dtset=dtset) 444 call rhotoxc(enxc,kxc,mpi_enreg,nfft,dtset%ngfft,& 445 & work,0,work,0,nkxc,nk3xc,non_magnetic_xc,n3xccc,option,dtset%paral_kgb,rhor,rprimd,strsxc,1,& 446 & vxc,vxcavg,xccc3d,xcdata,k3xc=k3xc,vhartr=vhartr) 447 ABI_DEALLOCATE(work) 448 449 ABI_DEALLOCATE(vhartr) 450 ABI_DEALLOCATE(vxc) 451 ABI_DEALLOCATE(xccc3d) 452 453 !Initialize finite difference calculation of the ddk 454 455 call status(0,dtfil%filstat,iexit,level,'call getshell ') 456 nkpt3 = 0 457 458 !Prepare first call to getkgrid (obtain number of k points in FBZ) 459 dum_kptrlatt(:,:) = dtset%kptrlatt(:,:) 460 dum_nshiftk = dtset%nshiftk 461 ABI_CHECK(dum_nshiftk<=210,"dum_nshiftk must be <= 210!") 462 dum_shiftk(:,:) = zero 463 dum_shiftk(:,1:dtset%nshiftk) = dtset%shiftk(:,1:dtset%nshiftk) 464 dum_vacuum(:) = 0 465 466 ABI_ALLOCATE(dum_kptns,(3,0)) 467 ABI_ALLOCATE(dum_wtk,(0)) 468 call getkgrid(0,0,dtset%iscf,dum_kptns,3,dum_kptrlatt,& 469 & rdum,dtset%nsym,0,nkpt3,dum_nshiftk,dtset%nsym,& 470 & rprimd,dum_shiftk,dtset%symafm,dtset%symrel,& 471 & dum_vacuum,dum_wtk) 472 ABI_DEALLOCATE(dum_kptns) 473 ABI_DEALLOCATE(dum_wtk) 474 475 write(std_out,*) 'nonlinear : nkpt, nkpt3 = ',nkpt,nkpt3 476 !call flush(6) 477 !jmb : malloc() problem with gcc461_openmpi under max2 : change order of allocations works ?!? 478 !allocate(kneigh(30,nkpt),kg_neigh(30,nkpt,3),mvwtk(30,nkpt)) 479 ABI_ALLOCATE(kg_neigh,(30,nkpt,3)) 480 ABI_ALLOCATE(mvwtk,(30,nkpt)) 481 ABI_ALLOCATE(kneigh,(30,nkpt)) 482 483 ABI_ALLOCATE(kptindex,(2,nkpt3)) 484 ABI_ALLOCATE(kpt3,(3,nkpt3)) 485 486 call getshell(gmet,kneigh,kg_neigh,kptindex,dtset%kptopt,& 487 & dtset%kptrlatt,dtset%kptns,kpt3,mkmem,mkmem_max,mpi_enreg,mvwtk,& 488 & nkpt,nkpt3,nneigh,dtset%nshiftk,rmet,rprimd,dtset%shiftk,dtset%wtk) 489 490 ABI_ALLOCATE(pwind,(mpw,nneigh,mkmem)) 491 ABI_ALLOCATE(cgindex,(nkpt,nsppol)) 492 ABI_ALLOCATE(mpi_enreg%kpt_loc2ibz_sp,(0:mpi_enreg%nproc-1,1:mkmem_max, 1:2)) 493 ABI_ALLOCATE(mpi_enreg%mkmem,(0:mpi_enreg%nproc-1)) 494 495 call status(0,dtfil%filstat,iexit,level,'call initmv ') 496 call initmv(cgindex,dtset,gmet,kg,kneigh,kg_neigh,kptindex,& 497 & kpt3,mband,mkmem,mpi_enreg,mpw,dtset%nband,nkpt,& 498 & nkpt3,nneigh,npwarr,nsppol,occ,pwind) 499 500 call status(0,dtfil%filstat,iexit,level,'call dfptnl_loop ') 501 call dfptnl_loop(blkflg,cg,cgindex,dtfil,dtset,d3lo,gmet,gprimd,gsqcut,& 502 & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,mband,mgfft,& 503 & mkmem,mkmem_max,dtset%mk1mem,mpert,mpi_enreg,mpw,mvwtk,natom,nfft,& 504 & nkpt,nkpt3,nkxc,nk3xc,nneigh,nspinor,nsppol,npwarr,occ,psps,pwind,& 505 & rfpert,rprimd,ucvol,xred) 506 507 write(message,'(a,a,a)')ch10,& 508 & ' --- Third order energy calculation completed --- ',ch10 509 call wrtout(ab_out,message,'COLL') 510 511 !Complete missing elements using symmetry operations 512 call status(0,dtfil%filstat,iexit,level,'call d3sym ') 513 call d3sym(blkflg,d3lo,indsym,mpert,natom,nsym,& 514 & symrec,dtset%symrel) 515 516 !Open the formatted derivative database file, and write the 517 !preliminary information 518 if (mpi_enreg%me == 0) then 519 call status(0,dtfil%filstat,iexit,level,'call ioddb8_ou') 520 521 dscrpt=' Note : temporary (transfer) database ' 522 523 call ddb_hdr_init(ddb_hdr,dtset,psps,pawtab,DDB_VERSION,dscrpt,& 524 & 1,xred=xred,occ=occ) 525 526 call ddb_hdr_open_write(ddb_hdr, dtfil%fnameabo_ddb, dtfil%unddb) 527 528 call ddb_hdr_free(ddb_hdr) 529 530 ! Call main output routine 531 call dfptnl_doutput(blkflg,d3lo,dtset%mband,mpert,dtset%nkpt,dtset%natom,dtset%ntypat,dtfil%unddb) 532 533 ! Close DDB 534 close(dtfil%unddb) 535 536 ! Compute tensors related to third-order derivatives 537 call nlopt(blkflg,carflg,d3lo,d3cart,gprimd,mpert,natom,rprimd,ucvol) 538 539 if ((d3e_pert1(natom+2)==1).and.(d3e_pert2(natom+2)==1).and. & 540 & (d3e_pert3(natom+2)==1)) then 541 542 flag = 1 543 i1pert = natom+2 544 545 d3cart(:,:,i1pert,:,i1pert,:,i1pert) = & 546 & d3cart(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 547 548 write(ab_out,*)ch10 549 write(ab_out,*)' Non-linear optical susceptibility tensor d (pm/V)' 550 write(ab_out,*)' in cartesian coordinates' 551 write(ab_out,*)' i1dir i2dir i3dir d' 552 553 do i1dir = 1, 3 554 do i2dir = 1, 3 555 do i3dir = 1, 3 556 write(ab_out,'(3(5x,i2),5x,f16.9)') i1dir,i2dir,i3dir,& 557 & d3cart(1,i1dir,i1pert,i2dir,i1pert,i3dir,i1pert) 558 if ((blkflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1).or.& 559 & (carflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1)) flag = 0 560 end do 561 end do 562 end do 563 564 if (flag == 0) then 565 write(message,'(a,a,a,a,a,a)')ch10,& 566 & ' dfptnl_doutput: WARNING -',ch10,& 567 & ' matrix of third-order energies incomplete,',ch10,& 568 & ' non-linear optical coefficients may be wrong.' 569 call wrtout(ab_out,message,'COLL') 570 call wrtout(std_out,message,'COLL') 571 end if 572 573 end if ! d3e_pert1,d3e_pert2,d3e_pert3 574 575 if (((maxval(d3e_pert1(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. & 576 & (d3e_pert3(natom+2)/=0)).or.& 577 ((maxval(d3e_pert2(1:natom))/=0).and.(d3e_pert1(natom+2)/=0).and. & 578 & (d3e_pert3(natom+2)/=0)).or.& 579 ((maxval(d3e_pert3(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. & 580 & (d3e_pert1(natom+2)/=0))) then 581 ! Perform a check if all relevant elements are available 582 583 flag = 1 584 do i1pert = 1, natom 585 do i1dir = 1, 3 586 do i2dir = 1, 3 587 do i3dir = 1, 3 588 if ((blkflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.& 589 (blkflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.& 590 (blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0 591 if ((carflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.& 592 (carflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.& 593 (carflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0 594 end do 595 end do 596 end do 597 end do 598 599 write(ab_out,*)ch10 600 write(ab_out,*)' First-order change in the electronic dielectric ' 601 write(ab_out,*)' susceptibility tensor (Bohr^-1)' 602 write(ab_out,*)' induced by an atomic displacement' 603 write(ab_out,*)' atom displacement' 604 605 do i1pert = 1,natom 606 do i1dir = 1,3 607 write(ab_out,'(1x,i4,9x,i2,3(3x,f16.9))')i1pert,i1dir,& 608 & d3cart(1,i1dir,i1pert,1,natom+2,:,natom+2) 609 write(ab_out,'(16x,3(3x,f16.9))')& 610 & d3cart(1,i1dir,i1pert,2,natom+2,:,natom+2) 611 write(ab_out,'(16x,3(3x,f16.9))')& 612 & d3cart(1,i1dir,i1pert,3,natom+2,:,natom+2) 613 end do 614 write(ab_out,*) 615 end do 616 617 if (flag == 0) then 618 write(message,'(a,a,a,a,a,a)')ch10,& 619 & ' dfptnl_doutput: WARNING -',ch10,& 620 & ' matrix of third-order energies incomplete,',ch10,& 621 & ' changes in the dielectric susceptibility may be wrong.' 622 call wrtout(ab_out,message,'COLL') 623 call wrtout(std_out,message,'COLL') 624 end if 625 626 end if ! d3e_pert1,d3e_pert2,d3e_pert3 627 end if ! mpi_enreg%me 628 629 ABI_DEALLOCATE(blkflg) 630 ABI_DEALLOCATE(carflg) 631 ABI_DEALLOCATE(cg) 632 ABI_DEALLOCATE(d3lo) 633 ABI_DEALLOCATE(d3cart) 634 ABI_DEALLOCATE(rhog) 635 ABI_DEALLOCATE(rhor) 636 ABI_DEALLOCATE(kneigh) 637 ABI_DEALLOCATE(kg_neigh) 638 ABI_DEALLOCATE(kptindex) 639 ABI_DEALLOCATE(kpt3) 640 ABI_DEALLOCATE(mvwtk) 641 ABI_DEALLOCATE(pwind) 642 ABI_DEALLOCATE(cgindex) 643 ABI_DEALLOCATE(d3e_pert1) 644 ABI_DEALLOCATE(d3e_pert2) 645 ABI_DEALLOCATE(d3e_pert3) 646 ABI_DEALLOCATE(rfpert) 647 ABI_DEALLOCATE(amass) 648 ABI_DEALLOCATE(grxc) 649 ABI_DEALLOCATE(kg) 650 ABI_DEALLOCATE(kxc) 651 ABI_DEALLOCATE(k3xc) 652 ABI_DEALLOCATE(indsym) 653 ABI_DEALLOCATE(npwarr) 654 ABI_DEALLOCATE(symrec) 655 ABI_DEALLOCATE(irrzon) 656 ABI_DEALLOCATE(phnons) 657 ABI_DEALLOCATE(mpi_enreg%kpt_loc2ibz_sp) 658 ABI_DEALLOCATE(mpi_enreg%mkmem) 659 660 !Clean the header 661 call hdr_free(hdr) 662 663 !As the etotal energy has no meaning here, we set it to zero 664 !(to avoid meaningless side-effects when comparing ouputs...) 665 etotal = zero 666 667 call status(0,dtfil%filstat,iexit,level,' exit ') 668 call timab(501,2,tsec) 669 670 end subroutine nonlinear