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