TABLE OF CONTENTS
- ABINIT/initmv
- ABINIT/m_nonlinear
- ABINIT/nonlinear
- m_nonlinear/dfptnl_doutput
- nonlinear/print_chi2
- nonlinear/print_dchidtau
ABINIT/initmv [ Functions ]
NAME
initmv
FUNCTION
Initialize finite difference calculation of the ddk im dfptnl_mv.f
INPUTS
dtset <type(dataset_type)> = all input variables in this dataset gmet(3,3) = reciprocal space metric tensor in bohr**-2 kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere kneigh(30,nkpt2) = index of the neighbours of each k-point kg_neigh(30,nkpt2,3) = necessary to construct the vector joining a k-point to its nearest neighbour in case of a single k-point, a line of k-points or a plane of k-points. See getshell.F90 for details kptindex(2,nkpt3)= index of the k-points in the reduced BZ related to a k-point in the full BZ kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ mband = maximum number of bands mkmem = number of k points which can fit in memory mpi_enreg = information about MPI parallelization mpw = maximum number of plane waves nband(nkpt*nsppol)=number of bands at each k point, for each polarization nkpt2 = number of k-points in the reduced BZ nkpt3 = number of k-points in the full BZ nneigh = total number of neighbours required to evaluate the finite difference formula npwarr(nkpt2)=number of planewaves at each k point nsppol = number of spin polarizations occ(mband*nkpt*nsppol) = occupation number for each band for each k
OUTPUT
cgindex(nkpt2,nsppol) = for each k-point, cgindex tores the location of the WF in the cg array
me = index of the current processor
ineigh = index of a neighbour
ikpt_loc = index of the iteration on ikpt on the current processor
ikpt_rbz = index of a k-point in the reduced BZ
pwind(mpw,nneigh,mkmem) = array used to compute the overlap matrix smat between k-points
(see initberry.f for more explanations)
SOURCE
1470 subroutine initmv(cgindex,dtset,gmet,kg,kneigh,kg_neigh,kptindex,& 1471 & kpt3,mband,mkmem,mpi_enreg,mpw,nband,nkpt2,& 1472 & nkpt3,nneigh,npwarr,nsppol,occ,pwind) 1473 1474 !Arguments ------------------------------------ 1475 !scalars 1476 integer,intent(in) :: mband,mkmem,mpw,nkpt2,nkpt3,nneigh,nsppol 1477 type(MPI_type),intent(inout) :: mpi_enreg 1478 type(dataset_type),intent(in) :: dtset 1479 !arrays 1480 integer,intent(in) :: kg(3,mpw*mkmem),kneigh(30,nkpt2),kg_neigh(30,nkpt2,3) 1481 integer,intent(in) :: nband(nkpt2*nsppol),npwarr(nkpt2),kptindex(2,nkpt3) 1482 integer,intent(out) :: cgindex(nkpt2,nsppol),pwind(mpw,nneigh,mkmem) 1483 real(dp),intent(in) :: gmet(3,3),kpt3(3,nkpt3),occ(mband*nkpt2*nsppol) 1484 1485 !Local variables------------------------------- 1486 !scalars 1487 integer :: flag,iband,icg,ierr,ikg,ikg1,ikpt,ikpt2,ikpt_loc,ikpt_rbz 1488 integer :: index,ineigh,ipw,isppol,jpw,nband_k,mband_occ,mband_occ_k,npw_k 1489 integer :: npw_k1,orig,spaceComm 1490 real(dp) :: ecut_eff,sdeg 1491 character(len=500) :: message 1492 !arrays 1493 integer :: dg(3) 1494 integer,allocatable :: kg1(:,:),kg1_k(:,:),npwar1(:),npwtot(:) 1495 real(dp) :: dk(3),dk_(3) 1496 real(dp),allocatable :: kpt1(:,:) 1497 1498 !************************************************************************ 1499 1500 !DEBUG 1501 !write(std_out,*)' initmv : enter ' 1502 !ENDDEBUG 1503 1504 if (xmpi_paral== 1) then 1505 ! BEGIN TF_CHANGES 1506 spaceComm=mpi_enreg%comm_cell 1507 ! END TF_CHANGES 1508 mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0 1509 mpi_enreg%mkmem(:) = 0 1510 end if 1511 1512 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 1513 ABI_MALLOC(kg1_k,(3,mpw)) 1514 ABI_MALLOC(kg1,(3,mkmem*mpw)) 1515 ABI_MALLOC(kpt1,(3,nkpt2)) 1516 ABI_MALLOC(npwar1,(nkpt2)) 1517 ABI_MALLOC(npwtot,(nkpt2)) 1518 kg1_k(:,:) = 0 1519 pwind(:,:,:) = 0 1520 cgindex(:,:) = 0 1521 1522 !Compute the number of occupied bands. 1523 !Check that it is the same for every k-point and that 1524 !nband(ikpt) is equal to this value 1525 1526 if (nsppol == 1) then 1527 sdeg = two 1528 else if (nsppol == 2) then 1529 sdeg = one 1530 end if 1531 1532 !DEBUG 1533 !write(std_out,*)' list of nband ' 1534 !do isppol = 1, nsppol 1535 !do ikpt = 1, nkpt2 1536 ! 1537 !nband_k = nband(ikpt + (isppol - 1)*nkpt2) 1538 !write(std_out,*)' isppol, ikpt, nband_k=',isppol, ikpt, nband_k 1539 !end do 1540 !end do 1541 !ENDDEBUG 1542 1543 index = 0 1544 do isppol = 1, nsppol 1545 do ikpt = 1, nkpt2 1546 1547 mband_occ_k = 0 1548 nband_k = nband(ikpt + (isppol - 1)*nkpt2) 1549 1550 do iband = 1, nband_k 1551 index = index + 1 1552 if (abs(occ(index) - sdeg) < tol8) mband_occ_k = mband_occ_k + 1 1553 end do 1554 1555 if (nband_k /= mband_occ_k) then 1556 write(message,'(a,a,a)')& 1557 & ' In a non-linear response calculation, nband must be equal ',ch10,& 1558 & ' to the number of valence bands.' 1559 ABI_ERROR(message) 1560 end if 1561 1562 ! Note that the number of bands can be different for spin up and spin down 1563 if (ikpt > 1) then 1564 if (mband_occ /= mband_occ_k) then 1565 message = 'The number of valence bands is not the same for every k-point' 1566 ABI_ERROR(message) 1567 end if 1568 else 1569 mband_occ = mband_occ_k 1570 end if 1571 1572 end do ! close loop over ikpt 1573 end do ! close loop over isppol 1574 1575 !Find the location of each wavefunction 1576 1577 icg = 0 1578 do isppol = 1, nsppol 1579 do ikpt = 1, nkpt2 1580 1581 1582 ! fab: inserted the shift due to the spin... 1583 1584 nband_k = dtset%nband(ikpt+(isppol - 1)*nkpt2) 1585 npw_k = npwarr(ikpt) 1586 1587 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,mpi_enreg%me)) cycle 1588 1589 cgindex(ikpt,isppol) = icg 1590 icg = icg + dtset%nspinor*npw_k*nband_k 1591 1592 end do 1593 end do 1594 1595 1596 !Build pwind 1597 1598 do ineigh = 1, nneigh 1599 1600 do ikpt = 1, nkpt2 1601 ikpt2 = kneigh(ineigh,ikpt) 1602 ikpt_rbz = kptindex(1,ikpt2) ! index of the k-point in the reduced BZ 1603 kpt1(:,ikpt) = dtset%kptns(:,ikpt_rbz) 1604 end do 1605 1606 ! Set up the basis sphere of plane waves at kpt1 1607 kg1(:,:) = 0 1608 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg1,& 1609 & kpt1,mkmem,dtset%nband,nkpt2,'PERS',mpi_enreg,mpw,& 1610 & npwar1,npwtot,dtset%nsppol) 1611 1612 ikg = 0 ; ikg1 = 0 ; ikpt_loc = 0 1613 1614 if(dtset%nsppol/=1)then 1615 if(mpi_enreg%nproc/=1)then 1616 message = ' At present, non-linear response calculations for spin-polarized system cannot be done in parallel.' 1617 ABI_ERROR(message) 1618 else 1619 isppol=1 1620 end if 1621 else 1622 isppol=1 1623 end if 1624 1625 do ikpt = 1, nkpt2 1626 1627 nband_k = dtset%nband(ikpt+(isppol - 1)*nkpt2) 1628 ikpt2 = kneigh(ineigh,ikpt) 1629 ikpt_rbz = kptindex(1,ikpt2) ! index of the k-point in the reduced BZ 1630 1631 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,-1,mpi_enreg%me)) cycle 1632 1633 ikpt_loc = ikpt_loc + 1 1634 1635 mpi_enreg%kpt_loc2ibz_sp(mpi_enreg%me, ikpt_loc, 1) = ikpt 1636 1637 flag = 0 1638 npw_k = npwarr(ikpt) 1639 npw_k1 = npwarr(ikpt_rbz) 1640 dk_(:) = kpt3(:,ikpt2) - dtset%kptns(:,ikpt) 1641 dk(:) = dk_(:) - nint(dk_(:)) + real(kg_neigh(ineigh,ikpt,:),dp) 1642 dg(:) = nint(dk(:) - dk_(:)) 1643 1644 1645 if (kptindex(2,ikpt2) == 0) then 1646 kg1_k(:,1:npw_k1) = kg1(:,ikg1+1:ikg1+npw_k1) 1647 if (dg(1)==0.and.dg(2)==0.and.dg(3)==0) flag = 1 1648 else 1649 kg1_k(:,1:npw_k1) = -1*kg1(:,ikg1+1:ikg1+npw_k1) 1650 end if 1651 1652 orig = 1 1653 do ipw = 1, npw_k 1654 do jpw = orig, npw_k1 1655 1656 if ((kg(1,ikg + ipw) == kg1_k(1,jpw) - dg(1)).and. & 1657 & (kg(2,ikg + ipw) == kg1_k(2,jpw) - dg(2)).and. & 1658 & (kg(3,ikg + ipw) == kg1_k(3,jpw) - dg(3))) then 1659 1660 pwind(ipw,ineigh,ikpt_loc) = jpw 1661 if (flag == 1) orig = jpw + 1 1662 exit 1663 1664 end if 1665 1666 end do 1667 end do 1668 1669 ikg = ikg + npw_k 1670 ikg1 = ikg1 + npw_k1 1671 1672 end do ! close loop over k-points 1673 1674 end do ! close loop over ineigh 1675 mpi_enreg%mkmem(mpi_enreg%me) = mkmem 1676 1677 call xmpi_sum(mpi_enreg%kpt_loc2ibz_sp,spaceComm,ierr) 1678 call xmpi_sum(mpi_enreg%mkmem,spaceComm,ierr) 1679 1680 !---------------------------------------------------------------------------- 1681 1682 ABI_FREE(kg1) 1683 ABI_FREE(kg1_k) 1684 ABI_FREE(kpt1) 1685 ABI_FREE(npwar1) 1686 ABI_FREE(npwtot) 1687 1688 end subroutine initmv
ABINIT/m_nonlinear [ Modules ]
NAME
m_nonlinear
FUNCTION
DFT calculations of non linear response functions.
COPYRIGHT
Copyright (C) 2002-2024 ABINIT group (MVeithen,MB,LB) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_nonlinear 23 24 use defs_basis 25 use defs_wvltypes 26 use m_wffile 27 use m_errors 28 use m_abicore 29 use m_xmpi 30 use m_hdr 31 use m_ebands 32 use m_xcdata 33 use m_dtset 34 use m_dtfil 35 36 use defs_datatypes, only : pseudopotential_type, ebands_t 37 use defs_abitypes, only : MPI_type 38 use m_fstrings, only : sjoin, itoa 39 use m_time, only : timab 40 use m_symtk, only : symmetrize_xred, littlegroup_q 41 use m_dynmat, only : d3sym, sytens 42 use m_ddb, only : ddb_type, nlopt 43 use m_ddb_hdr, only : ddb_hdr_type 44 use m_ioarr, only : read_rhor 45 use m_kg, only : getcut, kpgio, getph 46 use m_fft, only : fourdp 47 use m_kpts, only : getkgrid 48 use m_inwffil, only : inwffil 49 use m_spacepar, only : hartre, setsym 50 use m_pawfgr, only : pawfgr_type,pawfgr_init, pawfgr_destroy 51 use m_pawang, only : pawang_type, pawang_init, pawang_free 52 use m_pawrad, only : pawrad_type 53 use m_pawtab, only : pawtab_type,pawtab_get_lsize 54 use m_paw_an, only : paw_an_type, paw_an_init, paw_an_free, paw_an_nullify, paw_an_print 55 use m_paw_ij, only : paw_ij_type, paw_ij_init, paw_ij_free, paw_ij_nullify, paw_ij_print 56 use m_pawfgrtab, only : pawfgrtab_type, pawfgrtab_init, pawfgrtab_free 57 use m_pawrhoij, only : pawrhoij_type, pawrhoij_alloc, pawrhoij_free, pawrhoij_copy, & 58 & pawrhoij_bcast, pawrhoij_nullify, pawrhoij_inquire_dim 59 use m_pawdij, only : pawdij, symdij 60 use m_paw_finegrid,only : pawexpiqr 61 use m_pawxc, only : pawxc_get_nkxc 62 use m_paw_dmft, only : paw_dmft_type 63 use m_paw_sphharm, only : setsym_ylm 64 use m_paw_nhat, only : nhatgrid,pawmknhat 65 use m_paw_denpot, only : pawdenpot 66 use m_paw_init, only : pawinit,paw_gencond 67 use m_paw_tools, only : chkpawovlp 68 use m_mkrho, only : mkrho 69 use m_getshell, only : getshell 70 use m_pspini, only : pspini 71 use m_atm2fft, only : atm2fft 72 use m_rhotoxc, only : rhotoxc 73 use m_drivexc, only : check_kxc 74 use m_mpinfo, only : proc_distrb_cycle 75 use m_mklocl, only : mklocl 76 use m_common, only : setup1 77 use m_fourier_interpol, only : transgrid 78 use m_paw_occupancies, only : initrhoij 79 use m_paw_correlations, only : pawpuxinit 80 use m_mkcore, only : mkcore 81 use m_pead_nl_loop, only : pead_nl_loop 82 use m_dfptnl_loop, only : dfptnl_loop 83 84 implicit none 85 86 private
ABINIT/nonlinear [ Functions ]
NAME
nonlinear
FUNCTION
Primary routine for conducting DFT calculations of non linear response functions.
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) mpi_enreg=information about MPI pnarallelization 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
pawang <type(pawang_type)>=paw angular mesh and related data 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
NOTES
USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut)
for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ...
are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg)
for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ...
Total density, potentials, ...
are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used.
It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf)
are set equal to (nfft,ngfft,mgfft) in that case.
SOURCE
147 subroutine nonlinear(codvsn,dtfil,dtset,etotal,mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,xred) 148 149 !Arguments ------------------------------------ 150 !scalars 151 real(dp),intent(inout) :: etotal 152 character(len=8),intent(in) :: codvsn 153 type(MPI_type),intent(inout) :: mpi_enreg 154 type(datafiles_type),intent(in) :: dtfil 155 type(dataset_type),intent(inout) :: dtset 156 type(pawang_type),intent(inout) :: pawang 157 type(pseudopotential_type),intent(inout) :: psps 158 !arrays 159 integer,intent(out) :: npwtot(dtset%nkpt) 160 real(dp),intent(inout) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol),xred(3,dtset%natom) 161 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 162 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 163 164 !Local variables------------------------------- 165 !scalars 166 logical :: paral_atom,call_pawinit,qeq0 167 integer,parameter :: level=50,formeig=0,response=1,cplex1=1 168 integer :: ask_accurate,band_index,bantot,cplex,cplex_rhoij,dum_nshiftk,flag,gnt_option,gscase 169 integer :: has_dijnd,has_diju,has_kxc,has_k3xc 170 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 171 integer :: iatom,indx,iband,ider,idir,ierr,ifft,ikpt,ipert,isppol 172 integer :: ireadwf0,iscf_eff,ispden,itypat,izero,mcg,me,mgfftf,mkmem_max,mpert,my_natom 173 integer :: n1,n3xccc,natom,nband_k,nfftf,nfftot,nfftotf,nhatdim,nhatgrdim 174 integer :: nkpt_eff,nkpt_max,nkpt3,nkxc,nkxc1,nk3xc,nk3xc1,nneigh,ntypat,nsym1,nspden_rhoij,nzlmopt 175 integer :: optcut,optgr0,optgr1,optgr2,optrad,option,optorth 176 integer :: optatm,optdyfr,opteltfr,optgr,optstr,optv,optn,optn2 177 integer :: psp_gencond,pead,qphase_rhoij,rdwr,rdwrpaw,spaceworld,tim_mkrho,timrev 178 integer :: use_sym,usecprj,usexcnhat 179 logical :: is_dfpt=.true.,nmxc 180 real(dp),parameter :: k0(3)=(/zero,zero,zero/) 181 real(dp) :: boxcut,compch_fft,compch_sph,ecore,ecut_eff,ecutdg_eff,ecutf 182 real(dp) :: eei,epaw,epawdc,enxc,etot,fermie,fermih 183 real(dp) :: gsqcut,gsqcut_eff,gsqcutc_eff 184 real(dp) :: rdum,residm,ucvol,vxcavg 185 character(len=500) :: message 186 character(len=30) :: small_msg 187 character(len=fnlen) :: dscrpt 188 type(pawang_type) :: pawang1 189 type(ebands_t) :: bstruct 190 type(hdr_type) :: hdr,hdr_den 191 type(ddb_hdr_type) :: ddb_hdr 192 type(ddb_type) :: ddb 193 type(wffile_type) :: wffgs,wfftgs 194 type(wvl_data) :: wvl 195 type(xcdata_type) :: xcdata 196 !arrays 197 integer :: dum_kptrlatt(3,3),dum_vacuum(3),ngfft(18),ngfftf(18),perm(6),ii,theunit 198 integer,allocatable :: atindx(:),atindx1(:),blkflg(:,:,:,:,:,:),carflg(:,:,:,:,:,:),cgindex(:,:) 199 integer,allocatable :: flg_tmp(:,:,:,:,:,:) 200 integer,allocatable :: d3e_pert1(:),d3e_pert2(:),d3e_pert3(:) 201 integer,allocatable :: indsym(:,:,:),indsy1(:,:,:),irrzon(:,:,:),irrzon1(:,:,:) 202 integer,allocatable :: kg(:,:),kneigh(:,:),kg_neigh(:,:,:) 203 integer,allocatable :: kptindex(:,:),l_size_atm(:) 204 integer,allocatable :: npwarr(:),nattyp(:),pwind(:,:,:),rfpert(:,:,:,:,:,:) 205 integer,allocatable :: symq(:,:,:),symrec(:,:,:),symaf1(:),symrc1(:,:,:),symrl1(:,:,:) 206 real(dp) :: dum_gauss(0),dum_dyfrn(0),dum_dyfrv(0),dum_eltfrxc(0) 207 real(dp) :: dum_grn(0),dum_grv(0),dum_rhog(0),dum_vg(0) 208 real(dp) :: dum_shiftk(3,MAX_NSHIFTK),dummy6(6),other_dummy6(6),gmet(3,3),gprimd(3,3) 209 real(dp) :: qphon(3),rmet(3,3),rprimd(3,3),strsxc(6),tsec(2) 210 real(dp),allocatable :: cg(:,:),d3cart(:,:,:,:,:,:,:) 211 real(dp),allocatable :: d3etot(:,:,:,:,:,:,:),dum_kptns(:,:) 212 ! We need all these arrays instead of one because in Fortran the maximum number of dimensions is 7... 213 real(dp),allocatable :: d3e_1(:,:,:,:,:,:,:),d3cart_1(:,:,:,:,:,:,:) 214 real(dp),allocatable :: d3e_2(:,:,:,:,:,:,:),d3cart_2(:,:,:,:,:,:,:) 215 real(dp),allocatable :: d3e_3(:,:,:,:,:,:,:),d3cart_3(:,:,:,:,:,:,:) 216 real(dp),allocatable :: d3e_4(:,:,:,:,:,:,:),d3cart_4(:,:,:,:,:,:,:) 217 real(dp),allocatable :: d3e_5(:,:,:,:,:,:,:),d3cart_5(:,:,:,:,:,:,:) 218 real(dp),allocatable :: d3e_6(:,:,:,:,:,:,:),d3cart_6(:,:,:,:,:,:,:) 219 real(dp),allocatable :: d3e_7(:,:,:,:,:,:,:),d3cart_7(:,:,:,:,:,:,:) 220 real(dp),allocatable :: d3e_8(:,:,:,:,:,:,:),d3cart_8(:,:,:,:,:,:,:) 221 real(dp),allocatable :: d3e_9(:,:,:,:,:,:,:),d3cart_9(:,:,:,:,:,:,:) 222 real(dp),allocatable :: dum_wtk(:),dyfrlo_indx(:,:,:),dyfrx2(:,:,:),eigen0(:) 223 real(dp),allocatable :: grtn_indx(:,:),grxc(:,:),k3xc(:,:),kpt3(:,:),kxc(:,:) 224 real(dp),allocatable :: mvwtk(:,:),nhat(:,:),nhatgr(:,:,:),ph1d(:,:),ph1df(:,:),phnons(:,:,:),phnons1(:,:,:) 225 real(dp),allocatable :: rhog(:,:),rhor(:,:),rhowfg(:,:),rhowfr(:,:),tnons1(:,:) 226 real(dp),allocatable :: vhartr(:),vpsp(:),vtrial(:,:),vxc(:,:),work(:),xccc3d(:) 227 type(pawfgr_type) :: pawfgr 228 type(pawrhoij_type),allocatable :: pawrhoij(:),pawrhoij_read(:) 229 type(pawfgrtab_type),allocatable,save :: pawfgrtab(:) 230 type(paw_an_type),allocatable :: paw_an(:) 231 type(paw_ij_type),allocatable :: paw_ij(:) 232 type(paw_dmft_type) :: paw_dmft 233 234 ! *********************************************************************** 235 236 DBG_ENTER("COLL") 237 238 call timab(501,1,tsec) 239 240 !Structured debugging if dtset%prtvol==-level 241 if(dtset%prtvol==-level)then 242 write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' nonlinear : enter , debug mode ' 243 call wrtout(std_out,message,'COLL') 244 end if 245 246 !Check if the perturbations asked in the input file can be computed 247 248 if (((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert2_phon == 1)).or. & 249 & ((dtset%d3e_pert1_phon == 1).and.(dtset%d3e_pert3_phon == 1)).or. & 250 & ((dtset%d3e_pert2_phon == 1).and.(dtset%d3e_pert3_phon == 1))) then 251 write(message,'(7a)')& 252 & 'You have asked for a third-order derivative with respect to',ch10,& 253 & '2 or more atomic displacements.',ch10,& 254 & 'This is not allowed yet.',ch10,& 255 & 'Action : change d3e_pert1_phon, d3e_pert2_phon or d3e_pert3_phon in your input file.' 256 ABI_ERROR(message) 257 end if 258 259 !Computation of third order derivatives from PEAD (pead=1) or full DPFT formalism (pead=0): 260 pead = dtset%usepead 261 if (pead==0) then 262 write(message, '(2a)' ) ch10,'NONLINEAR : PEAD=0, full DFPT computation of third order derivatives' 263 call wrtout(ab_out,message,'COLL') 264 call wrtout(std_out,message,'COLL') 265 end if 266 267 !Some data for parallelism 268 nkpt_max=50;if(xmpi_paral==1)nkpt_max=-1 269 my_natom=mpi_enreg%my_natom 270 paral_atom=(my_natom/=dtset%natom) 271 if (paral_atom) then 272 ABI_BUG(" Nonlinear routine is not available yet with parallelization over atoms...") 273 end if 274 275 !Init spaceworld 276 spaceworld=mpi_enreg%comm_cell 277 me = xmpi_comm_rank(spaceworld) 278 279 !Define FFT grid(s) sizes (be careful !) 280 !See NOTES in the comments at the beginning of this file. 281 call pawfgr_init(pawfgr,dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf) 282 283 ntypat=psps%ntypat 284 natom=dtset%natom 285 nfftot=product(ngfft(1:3)) 286 nfftotf=product(ngfftf(1:3)) 287 288 !Define the set of admitted perturbations taking into account 289 !the possible permutations 290 mpert=natom+6 291 ABI_MALLOC(blkflg,(3,mpert,3,mpert,3,mpert)) 292 ABI_MALLOC(carflg,(3,mpert,3,mpert,3,mpert)) 293 ABI_MALLOC(rfpert,(3,mpert,3,mpert,3,mpert)) 294 ABI_MALLOC(d3e_pert1,(mpert)) 295 ABI_MALLOC(d3e_pert2,(mpert)) 296 ABI_MALLOC(d3e_pert3,(mpert)) 297 ABI_MALLOC(d3etot,(2,3,mpert,3,mpert,3,mpert)) 298 ABI_MALLOC(d3cart,(2,3,mpert,3,mpert,3,mpert)) 299 if (pead==0) then 300 ABI_MALLOC(d3e_1,(2,3,mpert,3,mpert,3,mpert)) 301 ABI_MALLOC(d3e_2,(2,3,mpert,3,mpert,3,mpert)) 302 ABI_MALLOC(d3e_3,(2,3,mpert,3,mpert,3,mpert)) 303 ABI_MALLOC(d3e_4,(2,3,mpert,3,mpert,3,mpert)) 304 ABI_MALLOC(d3e_5,(2,3,mpert,3,mpert,3,mpert)) 305 ABI_MALLOC(d3e_6,(2,3,mpert,3,mpert,3,mpert)) 306 ABI_MALLOC(d3e_7,(2,3,mpert,3,mpert,3,mpert)) 307 ABI_MALLOC(d3e_8,(2,3,mpert,3,mpert,3,mpert)) 308 ABI_MALLOC(d3e_9,(2,3,mpert,3,mpert,3,mpert)) 309 d3e_1(:,:,:,:,:,:,:) = 0_dp 310 d3e_2(:,:,:,:,:,:,:) = 0_dp 311 d3e_3(:,:,:,:,:,:,:) = 0_dp 312 d3e_4(:,:,:,:,:,:,:) = 0_dp 313 d3e_5(:,:,:,:,:,:,:) = 0_dp 314 d3e_6(:,:,:,:,:,:,:) = 0_dp 315 d3e_7(:,:,:,:,:,:,:) = 0_dp 316 d3e_8(:,:,:,:,:,:,:) = 0_dp 317 d3e_9(:,:,:,:,:,:,:) = 0_dp 318 if (dtset%nonlinear_info>0) then 319 ABI_MALLOC(flg_tmp,(3,mpert,3,mpert,3,mpert)) 320 ABI_MALLOC(d3cart_1,(2,3,mpert,3,mpert,3,mpert)) 321 ABI_MALLOC(d3cart_2,(2,3,mpert,3,mpert,3,mpert)) 322 ABI_MALLOC(d3cart_3,(2,3,mpert,3,mpert,3,mpert)) 323 ABI_MALLOC(d3cart_4,(2,3,mpert,3,mpert,3,mpert)) 324 ABI_MALLOC(d3cart_5,(2,3,mpert,3,mpert,3,mpert)) 325 ABI_MALLOC(d3cart_6,(2,3,mpert,3,mpert,3,mpert)) 326 ABI_MALLOC(d3cart_7,(2,3,mpert,3,mpert,3,mpert)) 327 ABI_MALLOC(d3cart_8,(2,3,mpert,3,mpert,3,mpert)) 328 ABI_MALLOC(d3cart_9,(2,3,mpert,3,mpert,3,mpert)) 329 end if 330 end if 331 blkflg(:,:,:,:,:,:) = 0 332 d3etot(:,:,:,:,:,:,:) = 0_dp 333 rfpert(:,:,:,:,:,:) = 0 334 d3e_pert1(:) = 0 ; d3e_pert2(:) = 0 ; d3e_pert3(:) = 0 335 336 if (dtset%d3e_pert1_phon==1) d3e_pert1(dtset%d3e_pert1_atpol(1):dtset%d3e_pert1_atpol(2))=1 337 if (dtset%d3e_pert2_phon==1) d3e_pert2(dtset%d3e_pert2_atpol(1):dtset%d3e_pert2_atpol(2))=1 338 if (dtset%d3e_pert3_phon==1) d3e_pert3(dtset%d3e_pert3_atpol(1):dtset%d3e_pert3_atpol(2))=1 339 if (dtset%d3e_pert1_elfd/=0) d3e_pert1(natom+2)=1 340 if (dtset%d3e_pert2_elfd/=0) d3e_pert2(natom+2)=1 341 if (dtset%d3e_pert3_elfd/=0) d3e_pert3(natom+2)=1 342 343 do i1pert = 1, mpert 344 do i1dir = 1, 3 345 do i2pert = 1, mpert 346 do i2dir = 1, 3 347 do i3pert = 1, mpert 348 do i3dir = 1, 3 349 perm(1) = & 350 & d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) & 351 & *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) & 352 & *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir) 353 perm(2) = & 354 & d3e_pert1(i1pert)*dtset%d3e_pert1_dir(i1dir) & 355 & *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) & 356 & *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir) 357 perm(3) = & 358 & d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) & 359 & *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) & 360 & *d3e_pert3(i3pert)*dtset%d3e_pert3_dir(i3dir) 361 perm(4) = & 362 & d3e_pert1(i2pert)*dtset%d3e_pert1_dir(i2dir) & 363 & *d3e_pert2(i3pert)*dtset%d3e_pert2_dir(i3dir) & 364 & *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir) 365 perm(5) = & 366 & d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) & 367 & *d3e_pert2(i2pert)*dtset%d3e_pert2_dir(i2dir) & 368 & *d3e_pert3(i1pert)*dtset%d3e_pert3_dir(i1dir) 369 perm(6) = & 370 & d3e_pert1(i3pert)*dtset%d3e_pert1_dir(i3dir) & 371 & *d3e_pert2(i1pert)*dtset%d3e_pert2_dir(i1dir) & 372 & *d3e_pert3(i2pert)*dtset%d3e_pert3_dir(i2dir) 373 if (sum(perm(:)) > 0) rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 374 end do 375 end do 376 end do 377 end do 378 end do 379 end do 380 381 ! call timab(134,2,tsec) 382 ! call timab(135,1,tsec) 383 384 !Do symmetry stuff 385 ABI_MALLOC(irrzon,(nfftot**(1-1/dtset%nsym),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 386 ABI_MALLOC(phnons,(2,nfftot**(1-1/dtset%nsym),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 387 ABI_MALLOC(indsym,(4,dtset%nsym,natom)) 388 ABI_MALLOC(symrec,(3,3,dtset%nsym)) 389 irrzon=0;indsym=0;symrec=0;phnons=zero 390 !If the density is to be computed by mkrho, need irrzon and phnons 391 iscf_eff=0;if(dtset%getden==0)iscf_eff=1 392 call setsym(indsym,irrzon,iscf_eff,natom,& 393 & nfftot,ngfft,dtset%nspden,dtset%nsppol,dtset%nsym,& 394 & phnons,dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred) 395 396 !Symmetrize atomic coordinates over space group elements: 397 call symmetrize_xred(natom,dtset%nsym,dtset%symrel,dtset%tnons,xred,indsym=indsym) 398 399 call sytens(indsym,mpert,natom,dtset%nsym,rfpert,symrec,dtset%symrel) 400 401 write(message, '(a,a,a,a,a)' ) ch10, & 402 & ' The list of irreducible elements of the Raman and non-linear',& 403 & ch10,' optical susceptibility tensors is:',ch10 404 call wrtout(ab_out,message,'COLL') 405 call wrtout(std_out,message,'COLL') 406 407 write(message,'(12x,a)')& 408 & 'i1pert i1dir i2pert i2dir i3pert i3dir' 409 call wrtout(ab_out,message,'COLL') 410 call wrtout(std_out,message,'COLL') 411 n1 = 0 412 do i1pert = 1, natom + 2 413 do i1dir = 1, 3 414 do i2pert = 1, natom + 2 415 do i2dir = 1,3 416 do i3pert = 1, natom + 2 417 do i3dir = 1, 3 418 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 419 n1 = n1 + 1 420 write(message,'(2x,i4,a,6(5x,i3))') n1,')', & 421 & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir 422 call wrtout(ab_out,message,'COLL') 423 call wrtout(std_out,message,'COLL') 424 else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-2) then 425 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 426 if (dtset%nonlinear_info>0) then 427 ! n1 = n1 + 1 428 write(message,'(2x,i4,a,6(5x,i3),a)') n1,')', & 429 & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir,' => must be zero, not computed' 430 call wrtout(ab_out,message,'COLL') 431 call wrtout(std_out,message,'COLL') 432 end if 433 else if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==-1) then 434 if (dtset%nonlinear_info>0) then 435 ! n1 = n1 + 1 436 write(message,'(2x,i4,a,6(5x,i3),a)') n1,')', & 437 & i1pert,i1dir,i2pert,i2dir,i3pert,i3dir,' => symmetric of an other element, not computed' 438 call wrtout(ab_out,message,'COLL') 439 call wrtout(std_out,message,'COLL') 440 end if 441 end if 442 end do 443 end do 444 end do 445 end do 446 end do 447 end do 448 write(message,'(a,a)') ch10,ch10 449 call wrtout(ab_out,message,'COLL') 450 call wrtout(std_out,message,'COLL') 451 452 ! For abipy : 453 if (dtset%paral_rf == -1) then 454 write(std_out,'(a)')"--- !IrredPerts" 455 write(std_out,'(a)')'# List of irreducible perturbations for nonlinear' 456 write(std_out,'(a)')'irred_perts:' 457 458 n1 = 0 459 do i1pert = 1, natom + 2 460 do i1dir = 1, 3 461 do i2pert = 1, natom + 2 462 do i2dir = 1, 3 463 do i3pert = 1, natom + 2 464 do i3dir = 1,3 465 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 466 n1 = n1 + 1 467 write(std_out,'(a,i0)')" - i1pert: ",i1pert 468 write(std_out,'(a,i0)')" i1dir: ",i1dir 469 write(std_out,'(a,i0)')" i2pert: ",i2pert 470 write(std_out,'(a,i0)')" i2dir: ",i2dir 471 write(std_out,'(a,i0)')" i3pert: ",i3pert 472 write(std_out,'(a,i0)')" i3dir: ",i3dir 473 end if 474 end do 475 end do 476 end do 477 end do 478 end do 479 end do 480 write(std_out,'(a)')"..." 481 ABI_ERROR_NODUMP("aborting now") 482 end if 483 484 !Set up for iterations 485 call setup1(dtset%acell_orig(1:3,1),bantot,dtset,& 486 ecutdg_eff,ecut_eff,gmet,gprimd,gsqcut_eff,gsqcutc_eff,& 487 ngfftf,ngfft,dtset%nkpt,dtset%nsppol,& 488 response,rmet,dtset%rprim_orig(1:3,1:3,1),rprimd,ucvol,psps%usepaw) 489 490 !Set up the basis sphere of planewaves 491 ABI_MALLOC(kg,(3,dtset%mpw*dtset%mkmem)) 492 ABI_MALLOC(npwarr,(dtset%nkpt)) 493 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg,& 494 & dtset%kptns,dtset%mkmem,dtset%nband,dtset%nkpt,'PERS',mpi_enreg,dtset%mpw,npwarr,npwtot,& 495 & dtset%nsppol) 496 497 !Recompute first large sphere cut-off gsqcut, without taking into account dilatmx 498 ecutf=dtset%ecut 499 if (psps%usepaw==1) then 500 ecutf=dtset%pawecutdg 501 call wrtout(std_out,ch10//' FFT (fine) grid used in SCF cycle:','COLL') 502 end if 503 call getcut(boxcut,ecutf,gmet,gsqcut,dtset%iboxcut,std_out,k0,ngfftf) 504 505 !Open and read pseudopotential files 506 ecore = 0_dp 507 call pspini(dtset,dtfil,ecore,psp_gencond,gsqcutc_eff,gsqcut_eff,pawrad,pawtab,& 508 & psps,rprimd,comm_mpi=mpi_enreg%comm_cell) 509 510 !Initialize band structure datatype 511 bstruct = ebands_from_dtset(dtset, npwarr) 512 513 !Initialize PAW atomic occupancies 514 if (psps%usepaw==1) then 515 ABI_MALLOC(pawrhoij,(my_natom)) 516 call pawrhoij_nullify(pawrhoij) 517 call initrhoij(dtset%pawcpxocc,dtset%lexexch,dtset%lpawu, & 518 & my_natom,natom,dtset%nspden,dtset%nspinor,dtset%nsppol,dtset%ntypat,& 519 & pawrhoij,dtset%pawspnorb,pawtab,cplex1,dtset%spinat,dtset%typat,& 520 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 521 else 522 ABI_MALLOC(pawrhoij,(0)) 523 end if 524 525 !Initialize header 526 gscase=0 527 call hdr_init(bstruct,codvsn,dtset,hdr,pawtab,gscase,psps,wvl%descr, & 528 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 529 530 !Update header, with evolving variables, when available 531 !Here, rprimd, xred and occ are available 532 etot=hdr%etot ; fermie=hdr%fermie ; fermih=hdr%fermih ; residm=hdr%residm 533 534 !If parallelism over atom, hdr is distributed 535 call hdr%update(bantot,etot,fermie,fermih,& 536 residm,rprimd,occ,pawrhoij,xred,dtset%amu_orig(:,1), & 537 comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab) 538 539 !Clean band structure datatype (should use it more in the future !) 540 call ebands_free(bstruct) 541 542 !Initialize wavefunction files and wavefunctions. 543 ireadwf0=1 544 545 mcg=dtset%mpw*dtset%nspinor*dtset%mband*dtset%mkmem*dtset%nsppol 546 ABI_MALLOC_OR_DIE(cg,(2,mcg), ierr) 547 548 ABI_MALLOC(eigen0,(dtset%mband*dtset%nkpt*dtset%nsppol)) 549 eigen0(:)=zero ; ask_accurate=1 550 optorth=0 551 552 call inwffil(ask_accurate,cg,dtset,dtset%ecut,ecut_eff,eigen0,dtset%exchn2n3d,& 553 & formeig,hdr,ireadwf0,dtset%istwfk,kg,dtset%kptns,& 554 & dtset%localrdwf,dtset%mband,mcg,dtset%mkmem,mpi_enreg,dtset%mpw,& 555 & dtset%nband,ngfft,dtset%nkpt,npwarr,dtset%nsppol,dtset%nsym,& 556 & occ,optorth,dtset%symafm,dtset%symrel,dtset%tnons,& 557 & dtfil%unkg,wffgs,wfftgs,dtfil%unwffgs,dtfil%fnamewffk,wvl) 558 559 !Close wffgs, if it was ever opened (in inwffil) 560 if (ireadwf0==1) then 561 call WffClose(wffgs,ierr) 562 end if 563 564 if (psps%usepaw==1.and.ireadwf0==1) then 565 ! if parallelism, pawrhoij is distributed, hdr%pawrhoij is not 566 call pawrhoij_copy(hdr%pawrhoij,pawrhoij,comm_atom=mpi_enreg%comm_atom,& 567 & mpi_atmtab=mpi_enreg%my_atmtab) 568 end if 569 570 ! call timab(135,2,tsec) 571 ! call timab(136,1,tsec) 572 573 !Report on eigen0 values ! Should use prteigrs.F90 574 write(message, '(a,a)' ) 575 call wrtout(std_out,ch10//' respfn : eigen0 array','COLL') 576 nkpt_eff=dtset%nkpt 577 if( (dtset%prtvol==0.or.dtset%prtvol==1.or.dtset%prtvol==2) .and. dtset%nkpt>nkpt_max ) nkpt_eff=nkpt_max 578 band_index=0 579 do isppol=1,dtset%nsppol 580 do ikpt=1,dtset%nkpt 581 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 582 if(ikpt<=nkpt_eff)then 583 write(message, '(a,i2,a,i5)' )' isppol=',isppol,', k point number',ikpt 584 call wrtout(std_out,message,'COLL') 585 do iband=1,nband_k,4 586 write(message, '(a,4es16.6)')' ',eigen0(iband+band_index:min(iband+3,nband_k)+band_index) 587 call wrtout(std_out,message,'COLL') 588 end do 589 else if(ikpt==nkpt_eff+1)then 590 write(message,'(a,a)' )' respfn : prtvol=0, 1 or 2, stop printing eigen0.',ch10 591 call wrtout(std_out,message,'COLL') 592 end if 593 band_index=band_index+nband_k 594 end do 595 end do 596 597 !Allocation for forces and atomic positions (should be taken away, also argument ... ) 598 ABI_MALLOC(grxc,(3,natom)) 599 600 !Examine the symmetries of the q wavevector 601 ABI_MALLOC(symq,(4,2,dtset%nsym)) 602 timrev=1 603 604 ! By default use symmetries. 605 use_sym = 1 606 if (dtset%prtgkk == 1)then 607 use_sym = 0 608 call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol,use_sym=use_sym) 609 else 610 call littlegroup_q(dtset%nsym,dtset%qptn,symq,symrec,dtset%symafm,timrev,prtvol=dtset%prtvol) 611 end if 612 613 614 615 !Generate an index table of atoms, in order for them to be used 616 !type after type. 617 ABI_MALLOC(atindx,(natom)) 618 ABI_MALLOC(atindx1,(natom)) 619 ABI_MALLOC(nattyp,(ntypat)) 620 indx=1 621 do itypat=1,ntypat 622 nattyp(itypat)=0 623 do iatom=1,natom 624 if(dtset%typat(iatom)==itypat)then 625 atindx(iatom)=indx 626 atindx1(indx)=iatom 627 indx=indx+1 628 nattyp(itypat)=nattyp(itypat)+1 629 end if 630 end do 631 end do 632 633 !Compute structure factor phases for current atomic pos: 634 ABI_MALLOC(ph1d,(2,3*(2*dtset%mgfft+1)*natom)) 635 ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*natom)) 636 call getph(atindx,natom,ngfft(1),ngfft(2),ngfft(3),ph1d,xred) 637 638 if (psps%usepaw==1.and.pawfgr%usefinegrid==1) then 639 call getph(atindx,natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred) 640 else 641 ph1df(:,:)=ph1d(:,:) 642 end if 643 644 qeq0=(dtset%qptn(1)**2+dtset%qptn(2)**2+dtset%qptn(3)**2<1.d-14) 645 if (.not.qeq0) then 646 ABI_BUG('NONLINEAR with dtset%qptn!=0 is not implemented yet') 647 end if 648 649 !PAW: 1- Initialize values for several arrays depending only on atomic data 650 !2- Check overlap 651 !3- Identify FFT points in spheres and compute g_l(r).Y_lm(r) (and exp(-i.q.r) if needed) 652 !4- Allocate PAW specific arrays 653 !5- Compute perturbed local potential inside spheres 654 !6- Eventually open temporary storage files 655 if(psps%usepaw==1) then 656 ! 1-Initialize values for several arrays depending only on atomic data 657 658 gnt_option=2 659 660 ! Test if we have to call pawinit 661 call paw_gencond(Dtset,gnt_option,"test",call_pawinit) 662 663 if (psp_gencond==1.or.call_pawinit) then 664 ! Some gen-cond have to be added... 665 call timab(553,1,tsec) 666 call pawinit(dtset%effmass_free,gnt_option,zero,zero,dtset%pawlcutd,dtset%pawlmix,& 667 & psps%mpsang,dtset%pawnphi,dtset%nsym,dtset%pawntheta,& 668 & pawang,pawrad,dtset%pawspnorb,pawtab,dtset%pawxcdev,dtset%ixc,dtset%usepotzero) 669 call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,& 670 & rprimd,symrec,pawang%zarot) 671 672 ! Update internal values 673 call paw_gencond(Dtset,gnt_option,"save",call_pawinit) 674 675 call timab(553,2,tsec) 676 else 677 if (pawtab(1)%has_kij ==1) pawtab(1:psps%ntypat)%has_kij =2 678 if (pawtab(1)%has_nabla==1) pawtab(1:psps%ntypat)%has_nabla=2 679 end if 680 psps%n1xccc=maxval(pawtab(1:psps%ntypat)%usetcore) 681 call setsym_ylm(gprimd,pawang%l_max-1,dtset%nsym,dtset%pawprtvol,rprimd,symrec,pawang%zarot) 682 call pawpuxinit(dtset%dmatpuopt,dtset%exchmix,dtset%f4of2_sla,dtset%f6of2_sla,& 683 & is_dfpt,dtset%jpawu,dtset%lexexch,dtset%lpawu,dtset%nspinor,ntypat,dtset%optdcmagpawu,pawang,dtset%pawprtvol,pawrad,& 684 & pawtab,dtset%upawu,dtset%usedmft,dtset%useexexch,dtset%usepawu) 685 compch_fft=-1.d5;compch_sph=-1.d5 686 usexcnhat=maxval(pawtab(:)%usexcnhat) 687 688 ! Note: many derivatives of cprj are needed and used a few times only, so for simplicity the 689 ! computation of all needed derivatives will be done on-the-fly. 690 usecprj=0 691 692 ! 2-Check overlap 693 call chkpawovlp(natom,psps%ntypat,dtset%pawovlp,pawtab,rmet,dtset%typat,xred) 694 ! 3-Identify FFT points in spheres and compute g_l(r).Y_lm(r) and exp(-i.q.r) 695 ABI_MALLOC(pawfgrtab,(my_natom)) 696 if (my_natom>0) then 697 call pawtab_get_lsize(pawtab,l_size_atm,my_natom,dtset%typat,& 698 & mpi_atmtab=mpi_enreg%my_atmtab) 699 call pawfgrtab_init(pawfgrtab,1,l_size_atm,pawrhoij(1)%nspden,dtset%typat,& 700 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 701 ABI_FREE(l_size_atm) 702 end if 703 optcut=0;optgr0=dtset%pawstgylm;optgr1=0;optgr2=0;optrad=1-dtset%pawstgylm 704 optgr1=dtset%pawstgylm 705 optgr2=dtset%pawstgylm 706 call nhatgrid(atindx1,gmet,my_natom,natom,nattyp,ngfftf,psps%ntypat,& 707 & optcut,optgr0,optgr1,optgr2,optrad,pawfgrtab,pawtab,rprimd,dtset%typat,ucvol,xred,& 708 & comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab ) 709 ABI_MALLOC(paw_an,(my_natom)) 710 ABI_MALLOC(paw_ij,(my_natom)) 711 call paw_an_nullify(paw_an) 712 call paw_ij_nullify(paw_ij) 713 has_kxc=0;nkxc1=0;cplex=1 714 has_dijnd=0;if(any(abs(dtset%nucdipmom)>tol8)) has_dijnd=1 715 has_diju=merge(0,1,dtset%usepawu==0) 716 has_kxc=1;nkxc1=2*dtset%nspden-1 ! LDA only 717 call pawxc_get_nkxc(nkxc1,dtset%nspden,dtset%xclevel) 718 has_k3xc=1; nk3xc1=3*min(dtset%nspden,2)-2 ! LDA only 719 call paw_an_init(paw_an,dtset%natom,dtset%ntypat,nkxc1,nk3xc1,dtset%nspden,cplex,dtset%pawxcdev,& 720 & dtset%typat,pawang,pawtab,has_vxc=1,has_vxc_ex=1,has_kxc=has_kxc,has_k3xc=has_k3xc,& 721 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 722 call paw_ij_init(paw_ij,cplex,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%pawspnorb,& 723 & natom,dtset%ntypat,dtset%typat,pawtab,has_dij=1,has_dijhartree=1,has_dijnd=has_dijnd,& 724 & has_dijso=1,has_dijU=has_diju,has_pawu_occ=1,has_exexch_pot=1,nucdipmom=dtset%nucdipmom,& 725 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 726 else ! PAW vs NCPP 727 usexcnhat=0;usecprj=0 728 ABI_MALLOC(paw_an,(0)) 729 ABI_MALLOC(paw_ij,(0)) 730 ABI_MALLOC(pawfgrtab,(0)) 731 end if 732 733 ABI_MALLOC(rhog,(2,nfftf)) 734 ABI_MALLOC(rhor,(nfftf,dtset%nspden)) 735 736 !Read ground-state charge density from diskfile in case getden /= 0 737 !or compute it from wfs that were read previously : rhor as well as rhog 738 739 if (dtset%getden /= 0 .or. dtset%irdden /= 0) then 740 ! Read rho1(r) from a disk file and broadcast data. 741 ! This part is not compatible with MPI-FFT (note single_proc=.True. below) 742 743 rdwr=1;rdwrpaw=psps%usepaw;if(ireadwf0/=0) rdwrpaw=0 744 if (rdwrpaw/=0) then 745 ABI_MALLOC(pawrhoij_read,(natom)) 746 call pawrhoij_nullify(pawrhoij_read) 747 call pawrhoij_inquire_dim(cplex_rhoij=cplex_rhoij,qphase_rhoij=qphase_rhoij,nspden_rhoij=nspden_rhoij,& 748 & nspden=dtset%nspden,spnorb=dtset%pawspnorb,cplex=cplex,cpxocc=dtset%pawcpxocc) 749 call pawrhoij_alloc(pawrhoij_read,cplex_rhoij,nspden_rhoij,dtset%nspinor,& 750 & dtset%nsppol,dtset%typat,qphase=qphase_rhoij,pawtab=pawtab) 751 else 752 ABI_MALLOC(pawrhoij_read,(0)) 753 end if 754 755 ! MT july 2013: Should we read rhoij from the density file ? 756 call read_rhor(dtfil%fildensin, cplex1, dtset%nspden, nfftf, ngfftf, rdwrpaw, mpi_enreg, rhor, & 757 hdr_den, pawrhoij_read, spaceworld, check_hdr=hdr) 758 call hdr_den%free() 759 760 if (rdwrpaw/=0) then 761 call pawrhoij_bcast(pawrhoij_read,hdr%pawrhoij,0,spaceworld) 762 call pawrhoij_free(pawrhoij_read) 763 end if 764 ABI_FREE(pawrhoij_read) 765 766 ! Compute up+down rho(G) by fft 767 ABI_MALLOC(work,(nfftf)) 768 work(:)=rhor(:,1) 769 call fourdp(1,rhog,work,-1,mpi_enreg,nfftf,1,ngfftf,0) 770 ABI_FREE(work) 771 772 else 773 izero=0 774 ! Obtain the charge density from read wfs 775 ! Be careful: in PAW, compensation density has to be added ! 776 tim_mkrho=4 777 paw_dmft%use_sc_dmft=0 ! respfn with dmft not implemented 778 paw_dmft%use_dmft=0 ! respfn with dmft not implemented 779 if (psps%usepaw==1) then 780 ABI_MALLOC(rhowfg,(2,dtset%nfft)) 781 ABI_MALLOC(rhowfr,(dtset%nfft,dtset%nspden)) 782 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 783 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhowfg,rhowfr,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 784 call transgrid(1,mpi_enreg,dtset%nspden,+1,1,1,dtset%paral_kgb,pawfgr,rhowfg,rhog,rhowfr,rhor) 785 ABI_FREE(rhowfg) 786 ABI_FREE(rhowfr) 787 else 788 call mkrho(cg,dtset,gprimd,irrzon,kg,mcg,& 789 & mpi_enreg,npwarr,occ,paw_dmft,phnons,rhog,rhor,rprimd,tim_mkrho,ucvol,wvl%den,wvl%wfs) 790 end if 791 end if ! getden 792 793 !In PAW, compensation density has eventually to be added 794 nhatgrdim=0;nhatdim=0 795 ABI_MALLOC(nhatgr,(0,0,0)) 796 if (psps%usepaw==1.and. ((usexcnhat==0).or.(dtset%getden==0).or.dtset%xclevel==2)) then 797 nhatdim=1 798 ABI_MALLOC(nhat,(nfftf,dtset%nspden)) 799 call timab(558,1,tsec) 800 nhatgrdim=0;if (dtset%xclevel==2.and.dtset%pawnhatxc>0) nhatgrdim=usexcnhat 801 ider=2*nhatgrdim 802 if (nhatgrdim>0) then 803 ABI_FREE(nhatgr) 804 ABI_MALLOC(nhatgr,(nfftf,dtset%nspden,3)) 805 end if 806 izero=0;cplex=1;ipert=0;idir=0;qphon(:)=zero 807 call pawmknhat(compch_fft,cplex,ider,idir,ipert,izero,gprimd,my_natom,natom,& 808 & nfftf,ngfftf,nhatgrdim,dtset%nspden,psps%ntypat,pawang,pawfgrtab,& 809 & nhatgr,nhat,pawrhoij,pawrhoij,pawtab,qphon,rprimd,ucvol,dtset%usewvl,xred, & 810 & mpi_atmtab=mpi_enreg%my_atmtab, comm_atom=mpi_enreg%comm_atom) 811 if (dtset%getden==0) then 812 rhor(:,:)=rhor(:,:)+nhat(:,:) 813 call fourdp(1,rhog,rhor(:,1),-1,mpi_enreg,nfftf,1,ngfftf,0) 814 end if 815 call timab(558,2,tsec) 816 else 817 ABI_MALLOC(nhat,(0,0)) 818 end if 819 820 !The GS irrzon and phnons were only needed to symmetrize the GS density 821 ABI_FREE(irrzon) 822 ABI_FREE(phnons) 823 824 !!jmb 2012 write(std_out,'(a)')' ' ! needed to make ibm6_xlf12 pass tests. No idea why this works. JWZ 5 Sept 2011 825 !!Will compute now the total potential 826 827 !Compute local ionic pseudopotential vpsp and core electron density xccc3d: 828 n3xccc=0;if (psps%n1xccc/=0) n3xccc=nfftf 829 ABI_MALLOC(xccc3d,(n3xccc)) 830 ABI_MALLOC(vpsp,(nfftf)) 831 832 eei = zero 833 if (psps%usepaw==1 .or. psps%nc_xccc_gspace==1) then 834 ! PAW or NC with nc_xccc_gspace: compute Vloc and core charge together in reciprocal space 835 call timab(562,1,tsec) 836 optatm=1;optdyfr=0;opteltfr=0;optgr=0;optstr=0;optv=1;optn=n3xccc/nfftf;optn2=1 837 call atm2fft(atindx1,xccc3d,vpsp,dum_dyfrn,dum_dyfrv,dum_eltfrxc,dum_gauss,gmet,gprimd,& 838 & dum_grn,dum_grv,gsqcut,mgfftf,psps%mqgrid_vl,natom,nattyp,nfftf,ngfftf,& 839 & ntypat,optatm,optdyfr,opteltfr,optgr,optn,optn2,optstr,optv,psps,pawtab,ph1df,psps%qgrid_vl,& 840 & dtset%qprtrb,dtset%rcut,dum_rhog,rprimd,dummy6,other_dummy6,ucvol,psps%usepaw,dum_vg,dum_vg,dum_vg,dtset%vprtrb,psps%vlspl) 841 call timab(562,2,tsec) 842 else 843 ! Norm-cons.: compute Vloc in reciprocal space and core charge in real space 844 option=1 845 ABI_MALLOC(dyfrlo_indx,(3,3,natom)) 846 ABI_MALLOC(grtn_indx,(3,natom)) 847 call mklocl(dtset,dyfrlo_indx,eei,gmet,gprimd,& 848 & grtn_indx,gsqcut,dummy6,mgfftf,mpi_enreg,natom,nattyp,& 849 & nfftf,ngfftf,dtset%nspden,ntypat,option,pawtab,ph1df,psps,& 850 & dtset%qprtrb,rhog,rhor,rprimd,ucvol,dtset%vprtrb,vpsp,wvl%descr,wvl%den,xred) 851 ABI_FREE(dyfrlo_indx) 852 ABI_FREE(grtn_indx) 853 if (psps%n1xccc/=0) then 854 ABI_MALLOC(dyfrx2,(3,3,natom)) 855 ABI_MALLOC(vxc,(0,0)) ! dummy 856 call mkcore(dummy6,dyfrx2,grxc,mpi_enreg,natom,nfftf,dtset%nspden,ntypat,& 857 & ngfftf(1),psps%n1xccc,ngfftf(2),ngfftf(3),option,rprimd,dtset%typat,ucvol,vxc,& 858 & psps%xcccrc,psps%xccc1d,xccc3d,xred) 859 ABI_FREE(dyfrx2) 860 ABI_FREE(vxc) ! dummy 861 end if 862 end if 863 864 !Set up hartree and xc potential. Compute kxc here. 865 ABI_MALLOC(vhartr,(nfftf)) 866 call hartre(1,gsqcut,dtset%icutcoul,psps%usepaw,mpi_enreg,nfftf,ngfftf,& 867 &dtset%nkpt,dtset%rcut,rhog,rprimd,dtset%vcutgeo,vhartr) 868 869 option=3 870 nkxc=2*dtset%nspden-1 ! LDA 871 if(dtset%xclevel==2.and.dtset%nspden==1) nkxc=7 ! non-polarized GGA 872 if(dtset%xclevel==2.and.dtset%nspden==2) nkxc=19 ! polarized GGA 873 nk3xc=3*dtset%nspden-2 874 call check_kxc(dtset%ixc,dtset%optdriver,check_k3xc=.true.) 875 ABI_MALLOC(kxc,(nfftf,nkxc)) 876 ABI_MALLOC(k3xc,(nfftf,nk3xc)) 877 ABI_MALLOC(vxc,(nfftf,dtset%nspden)) 878 879 call xcdata_init(xcdata,dtset=dtset) 880 nmxc=(dtset%usepaw==1.and.mod(abs(dtset%usepawu),10)==4) 881 call rhotoxc(enxc,kxc,mpi_enreg,nfftf,ngfftf,& 882 & nhat,nhatdim,nhatgr,nhatgrdim,nkxc,nk3xc,nmxc,n3xccc,option,rhor,& 883 & rprimd,strsxc,usexcnhat,vxc,vxcavg,xccc3d,xcdata,k3xc=k3xc,vhartr=vhartr) 884 885 !Compute local + Hxc potential, and subtract mean potential. 886 ABI_MALLOC(vtrial,(nfftf,dtset%nspden)) 887 do ispden=1,min(dtset%nspden,2) 888 do ifft=1,nfftf 889 vtrial(ifft,ispden)=vhartr(ifft)+vxc(ifft,ispden)+vpsp(ifft) 890 end do 891 end do 892 if (dtset%nspden==4) then 893 do ispden=3,4 894 do ifft=1,nfftf 895 vtrial(ifft,ispden)=vxc(ifft,ispden) 896 end do 897 end do 898 end if 899 ABI_FREE(vpsp) 900 ABI_FREE(vhartr) 901 902 if(dtset%prtvol==-level)then 903 call wrtout(std_out,' nonlinear : ground-state density and potential set up.','COLL') 904 end if 905 906 epaw = zero ; epawdc = zero 907 !PAW: compute Dij quantities (psp strengths) 908 if (psps%usepaw==1)then 909 cplex=1;ipert=0;option=1 910 nzlmopt=0;if (dtset%pawnzlm>0) nzlmopt=-1 911 call pawdenpot(compch_sph,epaw,epawdc,ipert,dtset%ixc,my_natom,natom,dtset%nspden,& 912 & ntypat,dtset%nucdipmom,nzlmopt,option,paw_an,paw_an,paw_ij,pawang,dtset%pawprtvol,& 913 & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,& 914 & dtset%spnorbscl,dtset%xclevel,dtset%xc_denpos,dtset%xc_taupos,ucvol,psps%znuclpsp, & 915 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 916 917 call timab(561,1,tsec) 918 call pawdij(cplex,dtset%enunit,gprimd,ipert,my_natom,natom,nfftf,nfftotf,& 919 & dtset%nspden,ntypat,paw_an,paw_ij,pawang,pawfgrtab,dtset%pawprtvol,& 920 & pawrad,pawrhoij,dtset%pawspnorb,pawtab,dtset%pawxcdev,k0,& 921 & dtset%spnorbscl,ucvol,dtset%cellcharge(1),vtrial,vxc,xred,nucdipmom=dtset%nucdipmom,& 922 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 923 call symdij(gprimd,indsym,ipert,my_natom,natom,dtset%nsym,ntypat,0,& 924 & paw_ij,pawang,dtset%pawprtvol,pawtab,rprimd,dtset%symafm,symrec,& 925 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 926 call timab(561,2,tsec) 927 end if 928 929 ABI_FREE(xccc3d) 930 931 ! Determine the subset of symmetry operations (nsym1 operations) 932 ! that leaves the perturbation invariant, and initialize corresponding arrays 933 ! symaf1, symrl1, tnons1 (and pawang1%zarot, if PAW).. 934 nsym1 = 1 935 ! symaf1_tmp(1) = 1 936 ! symrl1_tmp(:,:,1) = dtset%symrel(:,:,1) 937 ! tnons1_tmp(:,1) = 0_dp 938 ABI_MALLOC(indsy1,(4,nsym1,dtset%natom)) 939 ABI_MALLOC(symrc1,(3,3,nsym1)) 940 ABI_MALLOC(symaf1,(nsym1)) 941 ABI_MALLOC(symrl1,(3,3,nsym1)) 942 ABI_MALLOC(tnons1,(3,nsym1)) 943 symaf1(1)= 1 !symaf1_tmp(1:nsym1) 944 symrl1(:,:,1)= dtset%symrel(:,:,1) !symrl1_tmp(:,:,1:nsym1) 945 tnons1(:,1)= 0_dp !tnons1_tmp(:,1:nsym1) 946 ! ABI_FREE(symaf1_tmp) 947 ! ABI_FREE(symrl1_tmp) 948 ! ABI_FREE(tnons1_tmp) 949 950 ! Set up corresponding symmetry data 951 ABI_MALLOC(irrzon1,(dtset%nfft**(1-1/nsym1),2,(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 952 ABI_MALLOC(phnons1,(2,dtset%nfft**(1-1/nsym1),(dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4))) 953 call setsym(indsy1,irrzon1,1,dtset%natom,dtset%nfft,dtset%ngfft,dtset%nspden,dtset%nsppol,& 954 & nsym1,phnons1,symaf1,symrc1,symrl1,tnons1,dtset%typat,xred) 955 if (psps%usepaw==1) then 956 ! Allocate/initialize only zarot in pawang1 datastructure 957 call pawang_init(pawang1,0,0,pawang%l_max-1,0,0,nsym1,0,0,0,0) 958 call setsym_ylm(gprimd,pawang1%l_max-1,pawang1%nsym,0,rprimd,symrc1,pawang1%zarot) 959 end if 960 961 if (pead/=0) then 962 ! Initialize finite difference calculation of the ddk 963 964 nkpt3 = 0 965 966 ! Prepare first call to getkgrid (obtain number of k points in FBZ) 967 dum_kptrlatt(:,:) = dtset%kptrlatt(:,:) 968 dum_nshiftk = dtset%nshiftk 969 ABI_CHECK(dum_nshiftk <= MAX_NSHIFTK, sjoin("dum_nshiftk must be <= ", itoa(MAX_NSHIFTK))) 970 dum_shiftk(:,:) = zero 971 dum_shiftk(:,1:dtset%nshiftk) = dtset%shiftk(:,1:dtset%nshiftk) 972 dum_vacuum(:) = 0 973 974 ABI_MALLOC(dum_kptns,(3,0)) 975 ABI_MALLOC(dum_wtk,(0)) 976 call getkgrid(0,0,dtset%iscf,dum_kptns,3,dum_kptrlatt,& 977 & rdum,dtset%nsym,0,nkpt3,dum_nshiftk,dtset%nsym,& 978 & rprimd,dum_shiftk,dtset%symafm,dtset%symrel,& 979 & dum_vacuum,dum_wtk) 980 ABI_FREE(dum_kptns) 981 ABI_FREE(dum_wtk) 982 983 ! write(std_out,*) 'nonlinear : nkpt, nkpt3 = ',dtset%nkpt,nkpt3 984 !call flush(6) 985 !jmb : malloc() problem with gcc461_openmpi under max2 : change order of allocations works ?!? 986 !allocate(kneigh(30,nkpt),kg_neigh(30,nkpt,3),mvwtk(30,nkpt)) 987 ABI_MALLOC(kg_neigh,(30,dtset%nkpt,3)) 988 ABI_MALLOC(mvwtk,(30,dtset%nkpt)) 989 ABI_MALLOC(kneigh,(30,dtset%nkpt)) 990 991 ABI_MALLOC(kptindex,(2,nkpt3)) 992 ABI_MALLOC(kpt3,(3,nkpt3)) 993 994 call getshell(gmet,kneigh,kg_neigh,kptindex,dtset%kptopt,& 995 & dtset%kptrlatt,dtset%kptns,kpt3,dtset%mkmem,mkmem_max,mvwtk,& 996 & dtset%nkpt,nkpt3,nneigh,dtset%nshiftk,rmet,rprimd,dtset%shiftk,dtset%wtk, mpi_enreg%comm_cell) 997 998 ABI_MALLOC(pwind,(dtset%mpw,nneigh,dtset%mkmem)) 999 ABI_MALLOC(cgindex,(dtset%nkpt,dtset%nsppol)) 1000 ABI_MALLOC(mpi_enreg%kpt_loc2ibz_sp,(0:mpi_enreg%nproc-1,1:mkmem_max, 1:2)) 1001 ABI_MALLOC(mpi_enreg%mkmem,(0:mpi_enreg%nproc-1)) 1002 1003 call initmv(cgindex,dtset,gmet,kg,kneigh,kg_neigh,kptindex,& 1004 & kpt3,dtset%mband,dtset%mkmem,mpi_enreg,dtset%mpw,dtset%nband,dtset%nkpt,& 1005 & nkpt3,nneigh,npwarr,dtset%nsppol,occ,pwind) 1006 1007 call pead_nl_loop(blkflg,cg,cgindex,dtfil,dtset,d3etot,gmet,gprimd,gsqcut,& 1008 & hdr,kg,kneigh,kg_neigh,kptindex,kpt3,kxc,k3xc,dtset%mband,dtset%mgfft,& 1009 & dtset%mkmem,mkmem_max,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,mvwtk,natom,nfftf,& 1010 & dtset%nkpt,nkpt3,nkxc,nk3xc,nneigh,dtset%nspinor,dtset%nsppol,npwarr,occ,psps,pwind,& 1011 & rfpert,rprimd,ucvol,xred) 1012 1013 else ! pead=0 in this case 1014 1015 call dfptnl_loop(atindx,blkflg,cg,dtfil,dtset,d3etot,eigen0,gmet,gprimd,gsqcut,& 1016 & hdr,kg,kxc,k3xc,dtset%mband,dtset%mgfft,mgfftf,& 1017 & dtset%mkmem,dtset%mk1mem,mpert,mpi_enreg,dtset%mpw,natom,nattyp,ngfftf,nfftf,nhat,& 1018 & dtset%nkpt,nkxc,nk3xc,dtset%nspinor,dtset%nsppol,npwarr,occ,& 1019 & paw_an,paw_ij,pawang,pawang1,pawfgr,pawfgrtab,pawrad,pawrhoij,pawtab,& 1020 & ph1d,ph1df,psps,rfpert,rhog,rhor,rprimd,ucvol,usecprj,vtrial,vxc,xred,& 1021 & nsym1,indsy1,symaf1,symrc1,& 1022 & d3e_1,d3e_2,d3e_3,d3e_4,d3e_5,d3e_6,d3e_7,d3e_8,d3e_9) 1023 1024 !Complete missing elements using symmetry operations 1025 1026 if (dtset%nonlinear_info>0) then 1027 flg_tmp = blkflg 1028 call d3sym(flg_tmp,d3e_1,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1029 flg_tmp = blkflg 1030 call d3sym(flg_tmp,d3e_2,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1031 flg_tmp = blkflg 1032 call d3sym(flg_tmp,d3e_3,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1033 flg_tmp = blkflg 1034 call d3sym(flg_tmp,d3e_4,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1035 flg_tmp = blkflg 1036 call d3sym(flg_tmp,d3e_5,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1037 flg_tmp = blkflg 1038 call d3sym(flg_tmp,d3e_6,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1039 flg_tmp = blkflg 1040 call d3sym(flg_tmp,d3e_7,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1041 flg_tmp = blkflg 1042 call d3sym(flg_tmp,d3e_8,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1043 flg_tmp = blkflg 1044 call d3sym(flg_tmp,d3e_9,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1045 end if 1046 1047 end if ! end pead/=0 1048 1049 write(message,'(a,a,a)')ch10,& 1050 & ' --- Third order energy calculation completed --- ',ch10 1051 call wrtout(ab_out,message,'COLL') 1052 1053 1054 !Complete missing elements using symmetry operations 1055 call d3sym(blkflg,d3etot,indsym,mpert,natom,dtset%nsym,symrec,dtset%symrel) 1056 1057 if (mpi_enreg%me == 0) then 1058 1059 ! Write 3rd order derivatives in the output file 1060 call dfptnl_doutput(blkflg,d3etot,mpert) 1061 1062 ! Write the DDB file 1063 dscrpt=' Note : temporary (transfer) database ' 1064 call ddb_hdr%init(dtset,psps,pawtab,dscrpt,1,xred=xred,occ=occ) 1065 1066 call ddb%init(dtset, 1, mpert, with_d3E=.true.) 1067 1068 call ddb%set_d3matr(1, d3etot, blkflg) 1069 1070 call ddb%write_txt(ddb_hdr, dtfil%fnameabo_ddb) 1071 1072 call ddb_hdr%free() 1073 call ddb%free() 1074 1075 ! Compute tensors related to third-order derivatives 1076 call nlopt(blkflg,carflg,d3etot,d3cart,gprimd,mpert,natom,rprimd,ucvol) 1077 ! Note that the imaginary part is not transformed into cartesian coordinates 1078 if (pead==0.and.(dtset%nonlinear_info>0)) then 1079 call nlopt(blkflg,flg_tmp,d3e_1,d3cart_1,gprimd,mpert,natom,rprimd,ucvol) 1080 call nlopt(blkflg,flg_tmp,d3e_2,d3cart_2,gprimd,mpert,natom,rprimd,ucvol) 1081 call nlopt(blkflg,flg_tmp,d3e_3,d3cart_3,gprimd,mpert,natom,rprimd,ucvol) 1082 call nlopt(blkflg,flg_tmp,d3e_4,d3cart_4,gprimd,mpert,natom,rprimd,ucvol) 1083 call nlopt(blkflg,flg_tmp,d3e_5,d3cart_5,gprimd,mpert,natom,rprimd,ucvol) 1084 call nlopt(blkflg,flg_tmp,d3e_6,d3cart_6,gprimd,mpert,natom,rprimd,ucvol) 1085 call nlopt(blkflg,flg_tmp,d3e_7,d3cart_7,gprimd,mpert,natom,rprimd,ucvol) 1086 call nlopt(blkflg,flg_tmp,d3e_8,d3cart_8,gprimd,mpert,natom,rprimd,ucvol) 1087 call nlopt(blkflg,flg_tmp,d3e_9,d3cart_9,gprimd,mpert,natom,rprimd,ucvol) 1088 end if 1089 1090 if ((d3e_pert1(natom+2)==1).and.(d3e_pert2(natom+2)==1).and. & 1091 & (d3e_pert3(natom+2)==1)) then 1092 1093 flag = 1 1094 i1pert = natom+2 1095 1096 d3cart(:,:,i1pert,:,i1pert,:,i1pert) = & 1097 & d3cart(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 1098 1099 write(ab_out,*)ch10 1100 write(ab_out,*)' Non-linear optical susceptibility tensor d (pm/V)' 1101 write(ab_out,*)' in cartesian coordinates' 1102 write(ab_out,*)' i1dir i2dir i3dir d' 1103 1104 do i1dir = 1, 3 1105 do i2dir = 1, 3 1106 do i3dir = 1, 3 1107 write(ab_out,'(3(5x,i2),5x,f16.9)') i1dir,i2dir,i3dir,& 1108 & d3cart(1,i1dir,i1pert,i2dir,i1pert,i3dir,i1pert) 1109 if ((blkflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1).or.& 1110 & (carflg(i1dir,i1pert,i2dir,i1pert,i3dir,i1pert)/=1)) flag = 0 1111 end do 1112 end do 1113 end do 1114 1115 if (pead==0.and.(dtset%nonlinear_info>0)) then 1116 1117 d3cart_1(:,:,i1pert,:,i1pert,:,i1pert) = & 1118 & d3cart_1(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 1119 d3cart_2(:,:,i1pert,:,i1pert,:,i1pert) = & 1120 & d3cart_2(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 1121 d3cart_8(:,:,i1pert,:,i1pert,:,i1pert) = & 1122 & d3cart_8(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 1123 d3cart_9(:,:,i1pert,:,i1pert,:,i1pert) = & 1124 & d3cart_9(:,:,i1pert,:,i1pert,:,i1pert)*16*(pi**2)*(Bohr_Ang**2)*1.0d-8*eps0/e_Cb 1125 1126 theunit = ab_out 1127 1128 write(small_msg,'(a)') ' ** Total :' 1129 call print_chi2(d3cart,small_msg,theunit) 1130 1131 write(small_msg,'(a)') ' ** sum_psi1H1psi1 :' 1132 call print_chi2(d3cart_1,small_msg,theunit) 1133 1134 write(small_msg,'(a)') ' ** sum_lambda1psi1psi1 :' 1135 call print_chi2(d3cart_2,small_msg,theunit) 1136 1137 write(small_msg,'(a)') ' ** exc3 :' 1138 call print_chi2(d3cart_8,small_msg,theunit) 1139 1140 write(small_msg,'(a)') ' ** exc3_paw :' 1141 call print_chi2(d3cart_9,small_msg,theunit) 1142 1143 end if ! nonlinear_info > 0 1144 1145 if (flag == 0) then 1146 write(message,'(a,a,a,a,a,a)')ch10,& 1147 & ' dfptnl_doutput: WARNING -',ch10,& 1148 & ' matrix of third-order energies incomplete,',ch10,& 1149 & ' non-linear optical coefficients may be wrong, check input variables rfatpol and rfdir.' 1150 call wrtout(ab_out,message,'COLL') 1151 call wrtout(std_out,message,'COLL') 1152 end if 1153 1154 end if ! d3e_pert1,d3e_pert2,d3e_pert3 1155 1156 if (((maxval(d3e_pert1(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. & 1157 & (d3e_pert3(natom+2)/=0)).or.& 1158 ((maxval(d3e_pert2(1:natom))/=0).and.(d3e_pert1(natom+2)/=0).and. & 1159 & (d3e_pert3(natom+2)/=0)).or.& 1160 ((maxval(d3e_pert3(1:natom))/=0).and.(d3e_pert2(natom+2)/=0).and. & 1161 & (d3e_pert1(natom+2)/=0))) then 1162 ! Perform a check if all relevant elements are available 1163 1164 flag = 1 1165 do i1pert = 1, natom 1166 do i1dir = 1, 3 1167 do i2dir = 1, 3 1168 do i3dir = 1, 3 1169 if ((blkflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.& 1170 (blkflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.& 1171 (blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0 1172 if ((carflg(i1dir,i1pert,i2dir,natom+2,i3dir,natom+2) /= 1).or.& 1173 (carflg(i1dir,natom+2,i2dir,i1pert,i3dir,natom+2) /= 1).or.& 1174 (carflg(i1dir,natom+2,i2dir,natom+2,i3dir,i1pert) /= 1)) flag = 0 1175 end do 1176 end do 1177 end do 1178 end do 1179 1180 write(ab_out,*)ch10 1181 write(ab_out,*)' First-order change in the electronic dielectric ' 1182 write(ab_out,*)' susceptibility tensor (Bohr^-1)' 1183 write(ab_out,*)' induced by an atomic displacement' 1184 write(ab_out,*)' atom displacement' 1185 1186 do i1pert = 1,natom 1187 do i1dir = 1,3 1188 write(ab_out,'(1x,i4,9x,i2,3(3x,f16.9))')i1pert,i1dir,& 1189 & d3cart(1,i1dir,i1pert,1,natom+2,:,natom+2) 1190 write(ab_out,'(16x,3(3x,f16.9))')& 1191 & d3cart(1,i1dir,i1pert,2,natom+2,:,natom+2) 1192 write(ab_out,'(16x,3(3x,f16.9))')& 1193 & d3cart(1,i1dir,i1pert,3,natom+2,:,natom+2) 1194 end do 1195 write(ab_out,*) 1196 end do 1197 1198 if (flag == 0) then 1199 write(message,'(a,a,a,a,a,a)')ch10,& 1200 & ' dfptnl_doutput: WARNING -',ch10,& 1201 & ' matrix of third-order energies incomplete,',ch10,& 1202 & ' changes in the dielectric susceptibility may be wrong, check input variables rfatpol and rfdir.' 1203 call wrtout(ab_out,message,'COLL') 1204 call wrtout(std_out,message,'COLL') 1205 end if 1206 1207 if (pead==0.and.(dtset%nonlinear_info>0)) then 1208 theunit = ab_out 1209 1210 write(small_msg,'(a)') ' ** Total :' 1211 call print_dchidtau(d3cart,small_msg,theunit) 1212 1213 write(small_msg,'(a)') ' ** sum_psi1H1psi1 :' 1214 call print_dchidtau(d3cart_1,small_msg,theunit) 1215 1216 write(small_msg,'(a)') ' ** sum_lambda1psi1psi1 :' 1217 call print_dchidtau(d3cart_2,small_msg,theunit) 1218 1219 write(small_msg,'(a)') ' ** sum_lambda1psi0S1psi1 :' 1220 call print_dchidtau(d3cart_3,small_msg,theunit) 1221 1222 write(small_msg,'(a)') ' ** sum_psi0H2psi1a :' 1223 call print_dchidtau(d3cart_4,small_msg,theunit) 1224 1225 write(small_msg,'(a)') ' ** sum_psi0H2psi1b :' 1226 call print_dchidtau(d3cart_5,small_msg,theunit) 1227 1228 write(small_msg,'(a)') ' ** eHxc21_paw :' 1229 call print_dchidtau(d3cart_6,small_msg,theunit) 1230 1231 write(small_msg,'(a)') ' ** eHxc21_nhat :' 1232 call print_dchidtau(d3cart_7,small_msg,theunit) 1233 1234 write(small_msg,'(a)') ' ** exc3 :' 1235 call print_dchidtau(d3cart_8,small_msg,theunit) 1236 1237 write(small_msg,'(a)') ' ** exc3_paw :' 1238 call print_dchidtau(d3cart_9,small_msg,theunit) 1239 1240 end if ! nonlinear_info > 0 1241 1242 end if ! d3e_pert1,d3e_pert2,d3e_pert3 1243 end if ! mpi_enreg%me 1244 1245 ! TO OPTIMIZE DEALLOCATION ! 1246 if (pead/=0) then 1247 ABI_FREE(cgindex) 1248 ABI_FREE(kg_neigh) 1249 ABI_FREE(kneigh) 1250 ABI_FREE(kptindex) 1251 ABI_FREE(kpt3) 1252 ABI_FREE(mpi_enreg%kpt_loc2ibz_sp) 1253 ABI_FREE(mpi_enreg%mkmem) 1254 ABI_FREE(mvwtk) 1255 ABI_FREE(pwind) 1256 else 1257 if (dtset%nonlinear_info>0) then 1258 ABI_FREE(d3cart_1) 1259 ABI_FREE(d3cart_2) 1260 ABI_FREE(d3cart_3) 1261 ABI_FREE(d3cart_4) 1262 ABI_FREE(d3cart_5) 1263 ABI_FREE(d3cart_6) 1264 ABI_FREE(d3cart_7) 1265 ABI_FREE(d3cart_8) 1266 ABI_FREE(d3cart_9) 1267 ABI_FREE(flg_tmp) 1268 end if 1269 ABI_FREE(d3e_1) 1270 ABI_FREE(d3e_2) 1271 ABI_FREE(d3e_3) 1272 ABI_FREE(d3e_4) 1273 ABI_FREE(d3e_5) 1274 ABI_FREE(d3e_6) 1275 ABI_FREE(d3e_7) 1276 ABI_FREE(d3e_8) 1277 ABI_FREE(d3e_9) 1278 end if 1279 ABI_FREE(atindx) 1280 ABI_FREE(atindx1) 1281 ABI_FREE(blkflg) 1282 ABI_FREE(carflg) 1283 ABI_FREE(cg) 1284 ABI_FREE(d3cart) 1285 ABI_FREE(d3etot) 1286 ABI_FREE(d3e_pert1) 1287 ABI_FREE(d3e_pert2) 1288 ABI_FREE(d3e_pert3) 1289 ABI_FREE(eigen0) 1290 ABI_FREE(rhog) 1291 ABI_FREE(rhor) 1292 ABI_FREE(nhat) 1293 ABI_FREE(nhatgr) 1294 ABI_FREE(rfpert) 1295 ABI_FREE(grxc) 1296 ABI_FREE(kg) 1297 ABI_FREE(kxc) 1298 ABI_FREE(k3xc) 1299 ABI_FREE(indsym) 1300 ABI_FREE(indsy1) 1301 ABI_FREE(nattyp) 1302 ABI_FREE(npwarr) 1303 ABI_FREE(symrec) 1304 ABI_FREE(symrc1) 1305 ABI_FREE(symaf1) 1306 ABI_FREE(symrl1) 1307 ABI_FREE(tnons1) 1308 ABI_FREE(irrzon1) 1309 ABI_FREE(phnons1) 1310 ABI_FREE(symq) 1311 ABI_FREE(ph1d) 1312 ABI_FREE(ph1df) 1313 ABI_FREE(vtrial) 1314 ABI_FREE(vxc) 1315 call pawfgr_destroy(pawfgr) 1316 if (psps%usepaw==1) then 1317 call pawang_free(pawang1) 1318 call pawrhoij_free(pawrhoij) 1319 call paw_an_free(paw_an) 1320 call paw_ij_free(paw_ij) 1321 call pawfgrtab_free(pawfgrtab) 1322 end if 1323 ABI_FREE(pawrhoij) 1324 ABI_FREE(paw_an) 1325 ABI_FREE(paw_ij) 1326 ABI_FREE(pawfgrtab) 1327 1328 ! Clean the header 1329 call hdr%free() 1330 1331 !As the etotal energy has no meaning here, we set it to zero 1332 !(to avoid meaningless side-effects when comparing ouputs...) 1333 etotal = zero 1334 1335 call timab(501,2,tsec) 1336 1337 DBG_EXIT("COLL") 1338 1339 contains
m_nonlinear/dfptnl_doutput [ Functions ]
[ Top ] [ m_nonlinear ] [ Functions ]
NAME
dfptnl_doutput
FUNCTION
Write the matrix of third-order derivatives to the output file
INPUTS
blkflg(3,mpert,3,mpert,3,mpert)= ( 1 if the element of the 3dte has been calculated ; 0 otherwise ) d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE mpert =maximum number of ipert natom=Number of atoms ntypat=Number of type of atoms unddb = unit number for DDB output
NOTES
d3 holds the third-order derivatives before computing the permutations of the perturbations.
SOURCE
1715 subroutine dfptnl_doutput(blkflg,d3,mpert) 1716 1717 !Arguments ------------------------------- 1718 !scalars 1719 integer,intent(in) :: mpert 1720 !arrays 1721 integer,intent(in) :: blkflg(3,mpert,3,mpert,3,mpert) 1722 real(dp),intent(in) :: d3(2,3,mpert,3,mpert,3,mpert) 1723 1724 !Local variables ------------------------- 1725 !scalars 1726 integer :: i1dir,i1pert,i2dir,i2pert,i3dir,i3pert 1727 character(len=500) :: msg 1728 1729 !************************************************************************* 1730 1731 ! Write blok of third-order derivatives to ouput file 1732 1733 write(msg,'(a,a,a,a,a)')ch10,& 1734 ' Matrix of third-order derivatives (reduced coordinates)',ch10,& 1735 ' before computing the permutations of the perturbations',ch10 1736 call wrtout(ab_out,msg) 1737 1738 write(ab_out,*)' j1 j2 j3 matrix element' 1739 write(ab_out,*)' dir pert dir pert dir pert real part imaginary part' 1740 1741 do i1pert=1,mpert 1742 do i1dir=1,3 1743 do i2pert=1,mpert 1744 do i2dir=1,3 1745 do i3pert=1,mpert 1746 do i3dir=1,3 1747 1748 if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=0) then 1749 1750 write(ab_out,'(3(i4,i5),2f22.10)')& 1751 i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,& 1752 d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 1753 end if 1754 1755 end do 1756 end do 1757 end do 1758 end do 1759 end do 1760 end do 1761 1762 end subroutine dfptnl_doutput
nonlinear/print_chi2 [ Functions ]
[ Top ] [ nonlinear ] [ Functions ]
NAME
print_chi2
FUNCTION
Print a third derivative tensor. Used only in nonlinear
INPUTS
d3cart0 = the tensor to print msg = a short message printed before the tensor theunit = unit where the tensor is written
OUTPUT
SIDE EFFECTS
SOURCE
1360 subroutine print_chi2(d3cart0,msg,theunit) 1361 1362 integer,intent(in) :: theunit 1363 character(len=30) :: msg 1364 real(dp) :: elem1,elem2 1365 real(dp),intent(in) :: d3cart0(2,3,mpert,3,mpert,3,mpert) 1366 1367 ! ************************************************************************* 1368 1369 write(theunit,'(2a)') ch10,msg 1370 do i1dir = 1, 3 1371 do i2dir = 1, 3 1372 do i3dir = 1, 3 1373 elem1 = d3cart0(1,i1dir,natom+2,i2dir,natom+2,i3dir,natom+2) 1374 elem2 = d3cart0(2,i1dir,natom+2,i2dir,natom+2,i3dir,natom+2) 1375 write(theunit,'(3(5x,i2),5x,f16.9,2x,f16.9)') i1dir,i2dir,i3dir,elem1,elem2 1376 end do 1377 end do 1378 end do 1379 1380 end subroutine print_chi2
nonlinear/print_dchidtau [ Functions ]
[ Top ] [ nonlinear ] [ Functions ]
NAME
print_dchidtau
FUNCTION
Print a third derivative tensor. Used only in nonlinear
INPUTS
d3cart0 = the tensor to print msg = a short message printed before the tensor theunit = unit where the tensor is written
OUTPUT
SIDE EFFECTS
SOURCE
1401 subroutine print_dchidtau(d3cart0,msg,theunit) 1402 1403 integer,intent(in) :: theunit 1404 character(len=30) :: msg 1405 real(dp),intent(in) :: d3cart0(2,3,mpert,3,mpert,3,mpert) 1406 1407 ! ************************************************************************* 1408 1409 write(theunit,'(a)') msg 1410 do i1pert = 1,natom 1411 do i1dir = 1,3 1412 write(theunit,'(1x,i4,9x,i2,3(3x,f16.9),3(3x,f16.9))')i1pert,i1dir,& 1413 & d3cart0(1,i1dir,i1pert,1,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,1,natom+2,:,natom+2) 1414 write(theunit,'(16x,3(3x,f16.9),3(3x,f16.9))')& 1415 & d3cart0(1,i1dir,i1pert,2,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,2,natom+2,:,natom+2) 1416 write(theunit,'(16x,3(3x,f16.9),3(3x,f16.9))')& 1417 & d3cart0(1,i1dir,i1pert,3,natom+2,:,natom+2),d3cart0(2,i1dir,i1pert,3,natom+2,:,natom+2) 1418 end do 1419 end do 1420 1421 end subroutine print_dchidtau