Branch data Line data Source code
1 : : !!****m* ABINIT/m_psps
2 : : !! NAME
3 : : !! m_psps
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides method to allocate/free/initialize the pseudopotential_type object.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2014-2026 ABINIT group (XG,DC,MG)
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_psps
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : use m_xmpi
28 : : use m_nctk
29 : : use m_copy
30 : : use m_dtset
31 : : use netcdf
32 : :
33 : : use m_fstrings, only : itoa, sjoin, yesno, atoi
34 : : use m_io_tools, only : open_file
35 : : use m_matrix, only : matr3inv
36 : : use defs_datatypes, only : pspheader_type, pseudopotential_type, pseudopotential_gth_type, nctab_t
37 : : use m_paw_numeric, only : paw_spline
38 : : use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free, simp_gen
39 : : use m_pawpsp, only : pawpsp_cg
40 : : use m_parser, only : chkint_eq
41 : : use m_memeval, only : getdim_nloc, setmqgrid
42 : :
43 : : implicit none
44 : :
45 : : private
46 : :
47 : : ! Helper functions
48 : : public :: test_xml_xmlpaw_upf ! Test if a pseudo potential file is in XML, XML-PAW or in UPF format.
49 : :
50 : : !type, extends(pseudopotentials_base_t), public :: pseudopotentials_type
51 : : !contains
52 : : !end type pseudopotentials_base_type
53 : : !!***
54 : :
55 : : public :: psps_init_global ! Allocate and init all part of psps structure that are independent of a given dataset.
56 : : public :: psps_init_from_dtset ! Allocate and init all part of psps structure that are dependent of a given dataset.
57 : : public :: psps_free ! Deallocate all memory of psps structure.
58 : : public :: psps_copy ! Copy the psps structure.
59 : : public :: psps_print ! Print info on the pseudopotentials.
60 : : public :: psps_ncwrite_path ! Create a netcdf file and write psps data.
61 : : public :: psps_ncwrite ! Write psps data in an open netcdf file.
62 : : public :: psps_ncread ! Read psps data from an open netcdf file.
63 : :
64 : : !type, extends(nctab_base_t), public :: nctab_t
65 : : !contains
66 : : !end type nctab_t
67 : : !!***
68 : :
69 : : public :: nctab_init ! Create the object.
70 : : public :: nctab_free ! Free memory.
71 : : public :: nctab_copy ! Copy the object.
72 : : public :: nctab_eval_tvalespl ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space.
73 : : public :: nctab_eval_tvaletauspl ! Evaluate spline-fit of the atomic pseudo valence kinetic energy density in reciprocal space.
74 : : public :: nctab_eval_tcorespl ! Evalute spline-fit of the model core charge in reciprocal space.
75 : : public :: nctab_mixalch ! Mix the pseudopotential tables. Used for alchemical mixing.
76 : : !!***
77 : :
78 : : contains
79 : :
80 : : !!****f* m_psps/test_xml_xmlpaw_upf
81 : : !! NAME
82 : : !! test_xml_xmlpaw_upf
83 : : !!
84 : : !! FUNCTION
85 : : !! Test if a pseudo potential file is in XML, XML-PAW or in UPF format.
86 : : !!
87 : : !! INPUTS
88 : : !! path=Pseudopotential file
89 : : !!
90 : : !! OUTPUT
91 : : !! usexml=1 if XML file
92 : : !! xmlpaw=1 if PAW file in XML format
93 : : !! useupf=1 or 2 if UPF file.
94 : : !!
95 : : !! SOURCE
96 : :
97 : 2742 : subroutine test_xml_xmlpaw_upf(path, usexml, xmlpaw, useupf)
98 : :
99 : : !Arguments ------------------------------------
100 : : !scalars
101 : : character(len=*),intent(in) :: path
102 : : integer,intent(out) :: usexml, xmlpaw, useupf
103 : :
104 : : !Local variables-------------------------------
105 : : !scalars
106 : : integer :: temp_unit, ii
107 : : character(len=500) :: msg,errmsg
108 : : character(len=70) :: testxml
109 : : ! *************************************************************************
110 : :
111 : : ! Check if the file pseudopotential file is written in XML
112 : 2742 : usexml = 0; xmlpaw = 0; useupf = 0
113 : :
114 [ - + ]: 2742 : if (open_file(path,msg,newunit=temp_unit,form='formatted',status='old') /= 0) then
115 : 0 : ABI_ERROR(msg)
116 : : end if
117 [ + - ]: 2742 : rewind (unit=temp_unit,err=10,iomsg=errmsg)
118 : :
119 [ + - ]: 2742 : read(temp_unit, "(a)",err=10,iomsg=errmsg) testxml
120 [ + + ]: 2742 : if(testxml(1:5)=='<?xml')then
121 : 333 : usexml = 1
122 [ + - ]: 333 : read(temp_unit,*,err=10,iomsg=errmsg) testxml
123 [ + + ]: 333 : if(testxml(1:4)=='<paw') xmlpaw = 1
124 : : else
125 : 2409 : usexml = 0
126 [ + + ]: 2409 : if (testxml(1:4) == '<UPF') then
127 : : ! Make sure this is not UPF version >= 2
128 : : ! "<UPF version="2.0.1">
129 : 10 : ii = index(testxml, '"')
130 [ + - ]: 10 : if (ii /= 0) then
131 : 10 : useupf = atoi(testxml(ii+1:ii+1))
132 : : !if (useupf >= 2) then
133 : : ! ABI_ERROR(sjoin("UPF version >= 2 is not supported by Abinit. Use psp8 or psml format.", ch10, "Pseudo:", path))
134 : : !end if
135 : : else
136 : 0 : ABI_ERROR(sjoin("Cannot find version attributed in UPF file:", path))
137 : : end if
138 : :
139 : : end if
140 : : end if
141 : :
142 : : ! Check if pseudopotential file is a Q-espresso UPF1 file
143 [ + + ]: 2742 : if (useupf == 0) then
144 [ + - ]: 2732 : rewind (unit=temp_unit,err=10,iomsg=errmsg)
145 [ + - ]: 2732 : read(temp_unit,*,err=10,iomsg=errmsg) testxml ! just a string, no relation to xml.
146 [ + + ]: 2732 : if(testxml(1:9)=='<PP_INFO>')then
147 : 1 : useupf = 1
148 : : else
149 : 2731 : useupf = 0
150 : : end if
151 : : end if
152 : :
153 [ + - ]: 2742 : close(unit=temp_unit,err=10,iomsg=errmsg)
154 : :
155 : 2742 : return
156 : :
157 : : ! Handle IO error
158 : : 10 continue
159 : 0 : ABI_ERROR(errmsg)
160 : :
161 : : end subroutine test_xml_xmlpaw_upf
162 : : !!***
163 : :
164 : : !!****f* m_psps/psps_init_global
165 : : !! NAME
166 : : !! psps_init_global
167 : : !!
168 : : !! FUNCTION
169 : : !! Allocate and initialise all part of psps structure that are independent of a given dataset.
170 : : !!
171 : : !! INPUTS
172 : : !! npsp=the number of read pseudo files.
173 : : !! pspheads(npsp)=<type pspheader_type>all the important information from the
174 : : !! pseudopotential file header, as well as the psp file name
175 : : !!
176 : : !! SIDE EFFECTS
177 : : !! psps=<type pseudopotential_type>the pseudopotentials description
178 : : !!
179 : : !! SOURCE
180 : :
181 : 1410 : subroutine psps_init_global(psps, mtypalch, npsp, pspheads)
182 : :
183 : : !Arguments ------------------------------------
184 : : !scalars
185 : : class(pseudopotential_type),intent(inout) :: psps
186 : : integer,intent(in) :: mtypalch,npsp
187 : : !arrays
188 : : type(pspheader_type),intent(in) :: pspheads(npsp)
189 : :
190 : : !Local variables-------------------------------
191 : : integer :: ii, mpsang, n1xccc
192 : : ! *************************************************************************
193 : :
194 : : !Allocation of some arrays independent of the dataset
195 [ - + + - : 4230 : ABI_MALLOC(psps%filpsp,(npsp))
+ - + - -
+ - + -
+ ]
196 [ - + + - : 4230 : ABI_MALLOC(psps%pspcod,(npsp))
+ - + - -
+ - + -
+ ]
197 [ - + + - : 2820 : ABI_MALLOC(psps%pspdat,(npsp))
- + - + -
+ ]
198 [ - + + - : 2820 : ABI_MALLOC(psps%pspso,(npsp))
- + - + -
+ ]
199 [ - + + - : 2820 : ABI_MALLOC(psps%pspxc,(npsp))
- + - + -
+ ]
200 [ - + + - : 2820 : ABI_MALLOC(psps%title,(npsp))
- + - + -
+ ]
201 [ - + + - : 4230 : ABI_MALLOC(psps%zionpsp,(npsp))
+ - + - -
+ - + -
+ ]
202 [ - + + - : 2820 : ABI_MALLOC(psps%znuclpsp,(npsp))
- + - + -
+ ]
203 [ - + + - : 2820 : ABI_MALLOC(psps%epsatm,(npsp))
- + - + -
+ ]
204 : 1410 : call psp2params_init(psps%gth_params, npsp)
205 : :
206 [ + + ]: 3336 : psps%filpsp(1:npsp)=pspheads(1:npsp)%filpsp
207 [ + + ]: 3336 : psps%pspcod(1:npsp)=pspheads(1:npsp)%pspcod
208 [ + + ]: 3336 : psps%pspdat(1:npsp)=pspheads(1:npsp)%pspdat
209 [ + + ]: 3336 : psps%pspso(1:npsp)=pspheads(1:npsp)%pspso
210 [ + + ]: 3336 : psps%pspxc(1:npsp)=pspheads(1:npsp)%pspxc
211 [ + + ]: 3336 : psps%title(1:npsp)=pspheads(1:npsp)%title
212 [ + + ]: 3336 : psps%zionpsp(1:npsp)=pspheads(1:npsp)%zionpsp
213 [ + + ]: 3336 : psps%znuclpsp(1:npsp)=pspheads(1:npsp)%znuclpsp
214 : :
215 : : ! Transfer md5 checksum
216 [ - + + - : 4230 : ABI_MALLOC(psps%md5_pseudos, (npsp))
+ - + - -
+ - + -
+ ]
217 [ - + + + ]: 4746 : psps%md5_pseudos = pspheads(1:npsp)%md5_checksum
218 : : !Set values independant from dtset
219 : 1410 : psps%npsp = npsp
220 : : !Note that mpsang is the max of 1+lmax, with minimal value 1 (even for local psps, at present)
221 : 1410 : mpsang=1
222 : 1410 : n1xccc=pspheads(1)%xccc
223 [ + + ]: 3336 : do ii=1,psps%npsp
224 : 1926 : mpsang=max(pspheads(ii)%lmax+1,mpsang)
225 : 3336 : n1xccc=max(pspheads(ii)%xccc,n1xccc)
226 : : end do
227 : 1410 : psps%mpsang = mpsang
228 : 1410 : psps%n1xccc = n1xccc
229 : : ! Determine here whether the calculation is PAW
230 : : ! If paw, all pspcod necessarily are 7 or 17 (see iofn2)
231 : 1410 : psps%usepaw =0
232 [ + + ]: 1410 : if (pspheads(1)%pspcod==7.or.pspheads(1)%pspcod==17) psps%usepaw=1
233 : 1410 : psps%mtypalch = mtypalch
234 : :
235 : 1410 : end subroutine psps_init_global
236 : : !!***
237 : :
238 : : !----------------------------------------------------------------------
239 : :
240 : : !!****f* m_psps/psps_init_from_dtset
241 : : !! NAME
242 : : !! psps_init_from_dtset
243 : : !!
244 : : !! FUNCTION
245 : : !! Allocate and initialise all part of psps structure that are dependent of a given dataset.
246 : : !!
247 : : !! INPUTS
248 : : !! dtset=<type dataset_type>a given dataset
249 : : !! pspheads(npsp)=<type pspheader_type>all the important information from the
250 : : !! pseudopotential file header, as well as the psp file name
251 : : !!
252 : : !! SIDE EFFECTS
253 : : !! psps=<type pseudopotential_type>the pseudopotentials description
254 : : !!
255 : : !! SOURCE
256 : :
257 : 5340 : subroutine psps_init_from_dtset(psps, dtset, idtset, pspheads)
258 : :
259 : : !Arguments ------------------------------------
260 : : !scalars
261 : : class(pseudopotential_type),intent(inout) :: psps
262 : : integer,intent(in) :: idtset
263 : : type(dataset_type),intent(in) :: dtset
264 : : !arrays
265 : : type(pspheader_type),intent(in) :: pspheads(psps%npsp)
266 : :
267 : : !Local variables-------------------------------
268 : : !scalars
269 : : integer,save :: dimekb_old=-1,lmnmax_old=-1,lnmax_old=-1,mqgridff_old=0
270 : : integer,save :: mqgridvl_old=0,ntypat_old=-1,usepaw_old=-1
271 : : integer :: ipsp,lmnmax,lmnmaxso,lnmax,lnmaxso,newmqgrid,newmqgriddg,nptsgvec
272 : : integer :: changed,ii,itypat
273 : : real(dp) :: gprimd_orig(3,3)
274 : : ! *************************************************************************
275 : :
276 : 5340 : psps%optnlxccc = dtset%optnlxccc
277 : : !Determine the number of points needed in reciprocal space to represent the
278 : : !pseudopotentials (either set by hand from input variable or set automatically by abinit)
279 : 5340 : nptsgvec = 200 !This has to be chosen one and for all or else ??
280 : 5340 : newmqgrid = dtset%mqgrid
281 : 5340 : newmqgriddg = dtset%mqgriddg
282 : :
283 : : !JB:Which image to use ? I guess 1 always works
284 : 5340 : call matr3inv(dtset%rprimd_orig(:,:,1),gprimd_orig)
285 [ + - ]: 5340 : if ( dtset%usewvl == 0) then
286 : : call setmqgrid(newmqgrid,newmqgriddg,dtset%ecut*dtset%dilatmx**2,&
287 : 5340 : dtset%pawecutdg*dtset%dilatmx**2,gprimd_orig,nptsgvec,psps%usepaw)
288 : : else
289 : 0 : call setmqgrid(newmqgrid,newmqgriddg,one,one,gprimd_orig,nptsgvec,psps%usepaw)
290 : : end if
291 : 5340 : psps%mqgrid_ff = newmqgrid
292 [ + + ]: 5340 : if (psps%usepaw == 1) then
293 : 1213 : psps%mqgrid_vl = newmqgriddg
294 : : else
295 : 4127 : psps%mqgrid_vl = newmqgrid
296 : : end if
297 : :
298 : : !Determine the maximum number of projectors, for the set of pseudo atom
299 : : call getdim_nloc(lmnmax,lmnmaxso,lnmax,lnmaxso,dtset%mixalch_orig,dtset%nimage,psps%npsp,dtset%npspalch,&
300 : 5340 : dtset%ntypat,dtset%ntypalch,pspheads)
301 : :
302 : 5340 : psps%npspalch = dtset%npspalch
303 : 5340 : psps%ntypat = dtset%ntypat
304 : 5340 : psps%ntypalch = dtset%ntypalch
305 : 5340 : psps%ntyppure = dtset%ntyppure
306 : :
307 : : !Set the flag for reciprocal space or real space calculations
308 : 5340 : psps%vlspl_recipSpace = (dtset%icoulomb /= 1)
309 : 5340 : psps%positron = dtset%positron
310 : 5340 : psps%useylm = dtset%useylm
311 : 5340 : psps%usewvl = dtset%usewvl
312 : :
313 : : ! Define treatment of the model core density for NC pseudos.
314 : 5340 : psps%nc_xccc_gspace = dtset%nc_xccc_gspace
315 : :
316 [ + + ]: 5340 : if (idtset > 1) then
317 [ + - ]: 3930 : ABI_SFREE(psps%algalch)
318 [ + - ]: 3930 : ABI_SFREE(psps%mixalch)
319 : : end if
320 : :
321 [ + + + - : 16020 : ABI_MALLOC(psps%algalch,(psps%ntypalch))
+ - + + -
+ - + -
+ ]
322 [ + + + - : 21360 : ABI_MALLOC(psps%mixalch,(psps%npspalch,psps%ntypalch))
+ + + - +
- + + - +
- + - + ]
323 [ + + ]: 5360 : psps%algalch(1:psps%ntypalch)=dtset%algalch(1:psps%ntypalch)
324 : : !This value will be overwritten elsewhere in case there are different images ...
325 [ + + + + ]: 5400 : psps%mixalch(1:psps%npspalch,1:psps%ntypalch)=dtset%mixalch_orig(1:psps%npspalch,1:psps%ntypalch,1)
326 : :
327 : : !Set mpspso and psps%pspso
328 : : !Warning: mpspso might be different for each dataset.
329 : : ! mpspso not relevant in case of PAW.
330 : 5340 : psps%mpspso=1
331 [ + + ]: 12348 : do ipsp=1,dtset%npsp
332 [ + + ]: 12348 : if(dtset%nspinor==1)then
333 : 6477 : psps%pspso(ipsp)=0
334 : :
335 : : ! Ideally the following line should not exist, but at present, the space has to be booked
336 [ + + ]: 6477 : if(pspheads(ipsp)%pspso/=0)psps%mpspso=2
337 [ + + ]: 531 : else if (psps%usepaw==0) then
338 [ + + ]: 279 : if(dtset%so_psp(ipsp)/=1)then
339 : 129 : psps%pspso(ipsp)=dtset%so_psp(ipsp)
340 : : else
341 : 150 : psps%pspso(ipsp)=pspheads(ipsp)%pspso
342 : : end if
343 [ + + ]: 279 : if(psps%pspso(ipsp)/=0)psps%mpspso=2
344 [ + + ]: 279 : if(pspheads(ipsp)%pspso/=0)psps%mpspso=2
345 : : else
346 : 252 : psps%pspso(ipsp)=1+dtset%pawspnorb
347 : : end if
348 : : end do
349 : :
350 : : !Set mpssoang, lmnmax, lnmax
351 [ + + ]: 5340 : if(psps%mpspso==1)then
352 : 4350 : psps%mpssoang=psps%mpsang
353 : 4350 : psps%lmnmax =lmnmax
354 : 4350 : psps%lnmax =lnmax
355 : : else
356 : 990 : psps%mpssoang=2*psps%mpsang-1
357 : 990 : psps%lmnmax=lmnmaxso
358 : 990 : psps%lnmax=lnmaxso
359 : : end if
360 : :
361 : : !T. Rangel: for wvl + paw do not change psps%lmnmax
362 [ + + + - ]: 5340 : if (psps%useylm==0 .and. psps%usepaw/=1 ) then
363 : 3990 : psps%lmnmax=psps%lnmax
364 : : end if
365 : :
366 : : !Set dimekb
367 [ + + ]: 5340 : if (psps%usepaw==0) then
368 : 4127 : psps%dimekb=psps%lnmax
369 : : else
370 : 1213 : psps%dimekb=psps%lmnmax*(psps%lmnmax+1)/2
371 : : end if
372 : :
373 : : !The following arrays are often not deallocated before the end of the dtset loop
374 : : !and might keep their content from one dataset to the other, if the conditions are fulfilled
375 : 5340 : changed = 0
376 : :
377 [ + + + - : 5340 : if(dimekb_old/=psps%dimekb .or. ntypat_old/=dtset%ntypat .or. usepaw_old/=psps%usepaw) then
- + ]
378 : 1413 : changed = changed + 1
379 [ + + ]: 1413 : if(idtset/=1) then
380 [ + - ]: 3 : ABI_SFREE(psps%ekb)
381 : : end if
382 [ - + + - : 5652 : ABI_MALLOC(psps%ekb,(psps%dimekb,dtset%ntypat*(1-psps%usepaw)))
+ + + - +
- + + - +
- + - + ]
383 [ + + + + ]: 8277 : psps%ekb = zero
384 : 1413 : dimekb_old=psps%dimekb
385 : : end if
386 : :
387 [ + + - + ]: 5340 : if(lmnmax_old/=psps%lmnmax .or. ntypat_old/=dtset%ntypat)then
388 : 1457 : changed = changed + 1
389 [ + + ]: 1457 : if(idtset/=1) then
390 [ + - ]: 47 : ABI_SFREE(psps%indlmn)
391 : : end if
392 [ - + + - : 5828 : ABI_MALLOC(psps%indlmn,(6,psps%lmnmax,dtset%ntypat))
- + + - +
- + - - +
- + - + ]
393 [ + + + + : 107306 : psps%indlmn = zero
+ + ]
394 : 1457 : lmnmax_old=psps%lmnmax
395 : : end if
396 : :
397 [ + + + + : 5340 : if(mqgridff_old/=psps%mqgrid_ff .or. lnmax_old/=psps%lnmax .or. ntypat_old/=dtset%ntypat)then
- + ]
398 : 1415 : changed = changed + 1
399 [ + + ]: 1415 : if(idtset/=1) then
400 [ + - ]: 5 : ABI_SFREE(psps%ffspl)
401 [ + - ]: 5 : ABI_SFREE(psps%qgrid_ff)
402 : : end if
403 [ - + + - : 7075 : ABI_MALLOC(psps%ffspl,(psps%mqgrid_ff,2,psps%lnmax,dtset%ntypat))
- + + - -
+ + - + -
+ - - + -
+ - + ]
404 [ - + + - : 4245 : ABI_MALLOC(psps%qgrid_ff,(psps%mqgrid_ff))
+ - + - -
+ - + -
+ ]
405 [ + + + + : 51024372 : psps%ffspl = zero
+ + + + ]
406 [ + + ]: 4347019 : psps%qgrid_ff = zero
407 : 1415 : mqgridff_old=psps%mqgrid_ff
408 : 1415 : lnmax_old=psps%lnmax
409 : : end if
410 : :
411 [ + + - + ]: 5340 : if(mqgridvl_old/=psps%mqgrid_vl .or. ntypat_old/=dtset%ntypat)then
412 : 1413 : changed = changed + 1
413 [ + + ]: 1413 : if(idtset/=1) then
414 [ + - ]: 3 : ABI_SFREE(psps%qgrid_vl)
415 [ + - ]: 3 : ABI_SFREE(psps%vlspl)
416 [ + + ]: 3 : if (allocated(psps%nctab)) then
417 [ + + ]: 4 : do ii=1,size(psps%nctab)
418 : 4 : call nctab_free(psps%nctab(ii))
419 : : end do
420 [ + - + + : 4 : ABI_FREE(psps%nctab)
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
421 : : end if
422 : : end if
423 [ + + - + ]: 1413 : if (idtset/=1 .and. .not.psps%vlspl_recipSpace) then
424 [ # # ]: 0 : ABI_SFREE(psps%dvlspl)
425 : : end if
426 : :
427 [ - + + - : 4239 : ABI_MALLOC(psps%qgrid_vl,(psps%mqgrid_vl))
+ - + - -
+ - + -
+ ]
428 [ - + + - : 5652 : ABI_MALLOC(psps%vlspl,(psps%mqgrid_vl,2,dtset%ntypat))
- + + - +
- + - - +
- + - + ]
429 [ + + ]: 4362423 : psps%qgrid_vl = zero
430 [ + + + + : 11828058 : psps%vlspl = zero
+ + ]
431 : :
432 [ + + ]: 1413 : if (psps%usepaw == 0) then
433 : : ! If you change usepaw in the input, you will get what you deserve!
434 [ - + + - : 4375 : ABI_MALLOC(psps%nctab, (dtset%ntypat))
+ - + - -
+ - + - +
+ + ]
435 [ + + ]: 2347 : do itypat=1,dtset%ntypat
436 : 2347 : call nctab_init(psps%nctab(itypat), psps%mqgrid_vl, .False., .False.)
437 : : end do
438 : : end if
439 : :
440 [ - + ]: 1413 : if (.not.psps%vlspl_recipSpace) then
441 [ # # # # : 0 : ABI_MALLOC(psps%dvlspl,(psps%mqgrid_vl,2,dtset%ntypat))
# # # # #
# # # # #
# # # # ]
442 [ # # # # : 0 : psps%dvlspl = zero
# # ]
443 : : end if
444 : 1413 : mqgridvl_old=psps%mqgrid_vl
445 : : end if
446 : :
447 [ + + - + ]: 5340 : if(ntypat_old/=dtset%ntypat.or. usepaw_old/=psps%usepaw)then
448 : 1410 : changed = changed + 1
449 [ - + ]: 1410 : if(idtset/=1) then
450 [ # # ]: 0 : ABI_SFREE(psps%xccc1d)
451 [ # # ]: 0 : ABI_SFREE(psps%xcctau1d)
452 : : end if
453 [ + + + - : 5640 : ABI_MALLOC(psps%xccc1d,(psps%n1xccc*(1-psps%usepaw),6,dtset%ntypat))
- + + - +
- + + - +
- + - + ]
454 [ + + + - : 4230 : ABI_MALLOC(psps%xcctau1d,(psps%n1xccc*(1-psps%usepaw),6,dtset%ntypat))
- + + - -
+ - + -
+ ]
455 [ + + + + : 11914594 : psps%xccc1d = zero
+ + ]
456 [ + + + + : 11914594 : psps%xcctau1d = zero
+ + ]
457 : 1410 : usepaw_old=psps%usepaw
458 : : end if
459 : :
460 [ + + ]: 5340 : if(ntypat_old/=dtset%ntypat)then
461 : 1410 : changed = changed + 1
462 [ - + ]: 1410 : if(idtset/=1) then
463 [ # # ]: 0 : ABI_SFREE(psps%xcccrc)
464 [ # # ]: 0 : ABI_SFREE(psps%ziontypat)
465 [ # # ]: 0 : ABI_SFREE(psps%znucltypat)
466 : : end if
467 [ - + + - : 4230 : ABI_MALLOC(psps%xcccrc,(dtset%ntypat))
+ - + - -
+ - + -
+ ]
468 [ - + + - : 2820 : ABI_MALLOC(psps%znucltypat,(dtset%ntypat))
- + - + -
+ ]
469 [ - + + - : 2820 : ABI_MALLOC(psps%ziontypat,(dtset%ntypat))
- + - + -
+ ]
470 [ + + ]: 3328 : psps%xcccrc = zero
471 [ + + ]: 3328 : psps%znucltypat = zero
472 [ + + ]: 3328 : psps%ziontypat = zero
473 : 1410 : ntypat_old=dtset%ntypat
474 : : end if
475 : :
476 [ + + ]: 12328 : psps%ziontypat(:)=dtset%ziontypat(:)
477 : :
478 : 5340 : end subroutine psps_init_from_dtset
479 : : !!***
480 : :
481 : : !----------------------------------------------------------------------
482 : :
483 : : !!****f* m_psps/psps_free
484 : : !! NAME
485 : : !! psps_free
486 : : !!
487 : : !! FUNCTION
488 : : !! Deallocate all memory of psps structure.
489 : : !!
490 : : !! SOURCE
491 : :
492 : 7051 : subroutine psps_free(psps)
493 : :
494 : : !Arguments ------------------------------------
495 : : class(pseudopotential_type),intent(inout) :: psps
496 : :
497 : : !Local variables-------------------------------
498 : : integer :: ii
499 : : ! *************************************************************************
500 : :
501 : : !Allocation of some arrays independent of the dataset
502 [ + + ]: 7051 : ABI_SFREE(psps%filpsp)
503 [ + + ]: 7051 : ABI_SFREE(psps%pspcod)
504 [ + + ]: 7051 : ABI_SFREE(psps%pspdat)
505 [ + + ]: 7051 : ABI_SFREE(psps%pspso)
506 [ + + ]: 7051 : ABI_SFREE(psps%pspxc)
507 [ + + ]: 7051 : ABI_SFREE(psps%title)
508 [ + + ]: 7051 : ABI_SFREE(psps%algalch)
509 [ + + ]: 7051 : ABI_SFREE(psps%mixalch)
510 [ + + ]: 7051 : ABI_SFREE(psps%ekb)
511 [ + + ]: 7051 : ABI_SFREE(psps%indlmn)
512 [ + + ]: 7051 : ABI_SFREE(psps%ffspl)
513 [ + + ]: 7051 : ABI_SFREE(psps%qgrid_ff)
514 [ + + ]: 7051 : ABI_SFREE(psps%qgrid_vl)
515 [ + + ]: 7051 : ABI_SFREE(psps%vlspl)
516 [ - + ]: 7051 : ABI_SFREE(psps%dvlspl)
517 [ + + ]: 7051 : ABI_SFREE(psps%xccc1d)
518 [ + + ]: 7051 : ABI_SFREE(psps%xcctau1d)
519 [ + + ]: 7051 : ABI_SFREE(psps%xcccrc)
520 [ + + ]: 7051 : ABI_SFREE(psps%ziontypat)
521 [ + + ]: 7051 : ABI_SFREE(psps%zionpsp)
522 [ + + ]: 7051 : ABI_SFREE(psps%znucltypat)
523 [ + + ]: 7051 : ABI_SFREE(psps%znuclpsp)
524 [ + + ]: 7051 : ABI_SFREE(psps%md5_pseudos)
525 [ + + ]: 7051 : ABI_SFREE(psps%epsatm)
526 : :
527 : : ! Free types.
528 : 7051 : call psp2params_free(psps%gth_params)
529 : :
530 [ + + ]: 7051 : if (allocated(psps%nctab)) then
531 [ + + ]: 10972 : do ii=1,size(psps%nctab)
532 : 10972 : call nctab_free(psps%nctab(ii))
533 : : end do
534 [ + - + + : 10972 : ABI_FREE(psps%nctab)
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
535 : : end if
536 : :
537 : 7051 : end subroutine psps_free
538 : : !!***
539 : :
540 : : !----------------------------------------------------------------------
541 : :
542 : : !!****f* m_psps/psps_copy
543 : : !! NAME
544 : : !! psps_copy
545 : : !!
546 : : !! FUNCTION
547 : : !! Copy the psps structure.
548 : : !!
549 : : !! SOURCE
550 : :
551 : 3477 : subroutine psps_copy(pspsin, pspsout)
552 : :
553 : : !Arguments ------------------------------------
554 : : class(pseudopotential_type),intent(in) :: pspsin
555 : : class(pseudopotential_type),intent(inout) :: pspsout
556 : :
557 : : !Local variables-------------------------------
558 : : integer :: ii
559 : : ! *************************************************************************
560 : :
561 : : ! integer
562 : 3477 : pspsout%dimekb = pspsin%dimekb
563 : 3477 : pspsout%lmnmax = pspsin%lmnmax
564 : 3477 : pspsout%lnmax = pspsin%lnmax
565 : 3477 : pspsout%mproj = pspsin%mproj
566 : 3477 : pspsout%mpsang = pspsin%mpsang
567 : 3477 : pspsout%mpspso = pspsin%mpspso
568 : 3477 : pspsout%mpssoang = pspsin%mpssoang
569 : 3477 : pspsout%mqgrid_ff = pspsin%mqgrid_ff
570 : 3477 : pspsout%mqgrid_vl = pspsin%mqgrid_vl
571 : 3477 : pspsout%mtypalch = pspsin%mtypalch
572 : 3477 : pspsout%npsp = pspsin%npsp
573 : 3477 : pspsout%npspalch = pspsin%npspalch
574 : 3477 : pspsout%ntypat = pspsin%ntypat
575 : 3477 : pspsout%ntypalch = pspsin%ntypalch
576 : 3477 : pspsout%ntyppure = pspsin%ntyppure
577 : 3477 : pspsout%n1xccc = pspsin%n1xccc
578 : 3477 : pspsout%optnlxccc = pspsin%optnlxccc
579 : 3477 : pspsout%positron = pspsin%positron
580 : 3477 : pspsout%usepaw = pspsin%usepaw
581 : 3477 : pspsout%usewvl = pspsin%usewvl
582 : 3477 : pspsout%useylm = pspsin%useylm
583 : 3477 : pspsout%nc_xccc_gspace = pspsin%nc_xccc_gspace
584 : :
585 : : ! logical
586 : 3477 : pspsout%vlspl_recipSpace = pspsin%vlspl_recipSpace
587 : :
588 : : ! integer allocatable
589 [ + + - + ]: 3477 : if (allocated(pspsin%algalch)) call alloc_copy(pspsin%algalch, pspsout%algalch)
590 [ + - - + ]: 3477 : if (allocated(pspsin%indlmn)) call alloc_copy(pspsin%indlmn, pspsout%indlmn)
591 [ + - - + ]: 3477 : if (allocated(pspsin%pspdat)) call alloc_copy(pspsin%pspdat, pspsout%pspdat)
592 [ + - - + ]: 3477 : if (allocated(pspsin%pspcod)) call alloc_copy(pspsin%pspcod, pspsout%pspcod)
593 [ + - - + ]: 3477 : if (allocated(pspsin%pspso)) call alloc_copy(pspsin%pspso, pspsout%pspso)
594 [ + - - + ]: 3477 : if (allocated(pspsin%pspxc)) call alloc_copy(pspsin%pspxc, pspsout%pspxc)
595 : :
596 : : ! real allocatable
597 [ + - + + ]: 3477 : if (pspsin%dimekb > 0 .and. pspsin%usepaw==0) then
598 [ + - ]: 2799 : if (allocated(pspsin%ekb)) then
599 [ - + ]: 2799 : call alloc_copy( pspsin%ekb, pspsout%ekb)
600 : : end if
601 : : else
602 [ - + + - : 2712 : ABI_MALLOC(pspsout%ekb,(pspsout%dimekb,pspsout%ntypat * (1 - pspsout%usepaw)))
+ - - - +
- - + - +
- + - + ]
603 [ - + - - ]: 678 : pspsout%ekb = zero
604 : : end if
605 [ + + - + ]: 3477 : if (allocated(pspsin%ffspl)) call alloc_copy( pspsin%ffspl, pspsout%ffspl)
606 [ + + - + ]: 3477 : if (allocated(pspsin%mixalch)) call alloc_copy(pspsin%mixalch, pspsout%mixalch)
607 [ + + - + ]: 3477 : if (allocated(pspsin%qgrid_ff)) call alloc_copy(pspsin%qgrid_ff, pspsout%qgrid_ff)
608 [ + + - + ]: 3477 : if (allocated(pspsin%qgrid_vl)) call alloc_copy(pspsin%qgrid_vl, pspsout%qgrid_vl)
609 [ + + - + ]: 3477 : if (allocated(pspsin%vlspl)) call alloc_copy(pspsin%vlspl, pspsout%vlspl)
610 [ - + - - ]: 3477 : if (allocated(pspsin%dvlspl)) call alloc_copy(pspsin%dvlspl, pspsout%dvlspl)
611 : :
612 [ + - - + ]: 3477 : if (allocated(pspsin%ziontypat)) call alloc_copy(pspsin%ziontypat, pspsout%ziontypat)
613 [ + - - + ]: 3477 : if (allocated(pspsin%znucltypat)) call alloc_copy(pspsin%znucltypat, pspsout%znucltypat)
614 [ + + - + ]: 3477 : if (allocated(pspsin%epsatm)) call alloc_copy(pspsin%epsatm,pspsout%epsatm)
615 : :
616 : : ! GA: Could make a check on mtypalch here
617 [ + - - + ]: 3477 : if (allocated(pspsin%znuclpsp)) call alloc_copy(pspsin%znuclpsp, pspsout%znuclpsp)
618 [ + - - + ]: 3477 : if (allocated(pspsin%zionpsp)) call alloc_copy(pspsin%zionpsp, pspsout%zionpsp)
619 : :
620 [ + + ]: 3477 : if (pspsin%n1xccc > 0) then
621 [ + - - + ]: 2044 : if (allocated(pspsin%xcccrc)) call alloc_copy(pspsin%xcccrc, pspsout%xcccrc)
622 [ + - - + ]: 2044 : if (allocated(pspsin%xccc1d)) call alloc_copy(pspsin%xccc1d, pspsout%xccc1d)
623 [ + - - + ]: 2044 : if (allocated(pspsin%xcctau1d)) call alloc_copy(pspsin%xcctau1d, pspsout%xcctau1d)
624 : : end if
625 : :
626 : : ! allocate and copy character strings
627 [ - + + - : 10431 : ABI_MALLOC(pspsout%filpsp,(pspsout%npsp))
+ - + - -
+ - + -
+ ]
628 [ - + + - : 6954 : ABI_MALLOC(pspsout%title,(pspsout%npsp))
- + - + -
+ ]
629 [ - + + - : 10431 : ABI_MALLOC(pspsout%md5_pseudos,(pspsout%npsp))
+ - + - -
+ - + -
+ ]
630 [ + + ]: 7909 : do ii=1,pspsout%npsp
631 : 4432 : pspsout%filpsp(ii) = pspsin%filpsp(ii)
632 : 4432 : pspsout%title(ii) = pspsin%title(ii)
633 : 7909 : pspsout%md5_pseudos(ii) = pspsin%md5_pseudos(ii)
634 : : end do
635 : :
636 : : ! allocate and copy objects
637 [ + + ]: 3477 : if (allocated(pspsin%nctab)) then
638 [ - + + - : 11773 : ABI_MALLOC(pspsout%nctab,(pspsout%ntypat))
+ - + - -
+ - + - +
+ + ]
639 [ + + ]: 2800 : if (pspsin%usepaw==0) then
640 [ + + ]: 6170 : do ii=1,pspsout%ntypat
641 : 6170 : call nctab_copy(pspsin%nctab(ii), pspsout%nctab(ii))
642 : : end do
643 : : end if
644 : : end if
645 : :
646 : 3477 : call psp2params_copy(pspsin%gth_params, pspsout%gth_params)
647 : :
648 : 3477 : end subroutine psps_copy
649 : : !!***
650 : :
651 : : !----------------------------------------------------------------------
652 : :
653 : : !!****f* m_psps/psps_print
654 : : !! NAME
655 : : !! psps_print
656 : : !!
657 : : !! FUNCTION
658 : : !! Print the content of a pseudopotential_type derived type
659 : : !!
660 : : !! INPUTS
661 : : !! psps=<type pseudopotential_type>=Info on the pseudopotentials.
662 : : !! units=unit numbers for output
663 : : !! prtvol(optional)=verbosity level
664 : : !! mode_paral(optional): either "COLL" or "PERS"
665 : : !!
666 : : !! OUTPUT
667 : : !! Only writing
668 : : !!
669 : : !! SOURCE
670 : :
671 [ - + ]: 2287 : subroutine psps_print(psps, units, prtvol, mode_paral)
672 : :
673 : : !Arguments ------------------------------------
674 : : !scalars
675 : : class(pseudopotential_type),intent(in) :: psps
676 : : integer,intent(in) :: units(:)
677 : : integer,intent(in),optional :: prtvol
678 : : character(len=4),intent(in),optional :: mode_paral
679 : :
680 : : !Local variables-------------------------------
681 : : !scalars
682 : : integer :: ips,ipsp_alch,ityp_alch,itypat,my_prtvol
683 : : character(len=4) :: mode
684 : : character(len=500) :: msg
685 : : !arrays
686 : : ! *************************************************************************
687 : :
688 : : ! Provide defaults
689 [ + - ]: 2287 : my_prtvol=0; if (present(prtvol)) my_prtvol=prtvol
690 [ - + ]: 2287 : mode='COLL'; if (present(mode_paral)) mode=mode_paral
691 : :
692 : : ! General info including spin-orbit
693 : 2287 : call wrtout(units,' ==== Info on pseudopotentials ==== ', mode)
694 : :
695 : 4033 : select case (psps%usepaw)
696 : : case (0)
697 : 1746 : call wrtout(units,' Norm-conserving pseudopotentials ', mode)
698 : : !call wrtout(units, sjoin(' Max number of Kleinman-Bylander energies ', itoa(psps%dimekb)), mode)
699 : : !do itypat=1,psps%ntypat
700 : : ! write(msg,'(a,i4,a,f9.4)')' Type ',itypat,' K-B energies ',(psps%ekb(ikbe,itypat),ikbe=1,psps%dimekb)
701 : : !end do
702 : : case (1)
703 : 541 : write(msg,'(a)')
704 : 541 : call wrtout(units,' PAW calculation', mode)
705 : : !call wrtout(units,sjoin(' Max number of D_ij coefficients ', itoa(psps%dimekb)), mode)
706 : : case default
707 [ + + - ]: 2287 : ABI_ERROR(sjoin("Invalid usepaw: ", itoa(psps%usepaw)))
708 : : end select
709 : :
710 : : !select case (psps%positron)
711 : : !case (0)
712 : : ! call wrtout(units, ' Standard Electron Calculation ', mode)
713 : : !case (1,2)
714 : : ! write(msg,'(a,i0)')' Positron Calculation with positron .. ',psps%positron
715 : : ! call wrtout(units,msg,mode)
716 : : !case default
717 : : ! ABI_ERROR(sjoin("Invalid positron: ", itoa(psps%positron)))
718 : : !end select
719 : :
720 : : write(msg,'(a,i4,2a,i4)')&
721 : 2287 : ' Number of pseudopotentials .. ',psps%npsp,ch10,&
722 : 4574 : ' Number of types of atoms .. ',psps%ntypat
723 : 2287 : call wrtout(units,msg,mode)
724 : :
725 [ + + ]: 2287 : if (psps%usepaw==0) then
726 : 3159 : select case (psps%mpspso)
727 : : case (1)
728 : 1413 : call wrtout(units,' Scalar calculation (no spin-orbit term) ',mode)
729 : : case (2)
730 : : write(msg,'(3a,i3)')&
731 : 333 : ' Calculation with spin-orbit coupling ',ch10,&
732 : 666 : ' Max number of channels (spin-orbit included) ',psps%mpssoang
733 : 333 : call wrtout(units,msg,mode)
734 [ + + ]: 786 : do itypat=1,psps%ntypat
735 [ + - ]: 786 : if (psps%pspso(itypat) /= 1) then
736 : : write(msg,'(a,i4,a,i2,a)')&
737 : 453 : ' - Atom type ',itypat,' has spin-orbit characteristics (pspso= ',psps%pspso(itypat),")"
738 : 453 : call wrtout(units,msg,mode)
739 : : end if
740 : : end do
741 : : case default
742 [ + + - ]: 1746 : ABI_ERROR(sjoin("Invalid mpspso: ", itoa(psps%mpspso)))
743 : : end select
744 : : else
745 [ + + ]: 1312 : select case (maxval(psps%pspso))
746 : : case (0,1)
747 : 481 : msg=' Scalar calculation (no spin-orbit term) '
748 : : case (2)
749 [ + + - ]: 541 : msg=' Calculation with spin-orbit coupling '
750 : : end select
751 : 541 : call wrtout(units,msg,mode)
752 : : end if
753 : :
754 : : ! Info on nonlocal part
755 : 3974 : select case (psps%useylm)
756 : : case (0)
757 : 1687 : msg = ' Nonlocal part applied using Legendre polynomials '
758 : : case (1)
759 : 600 : msg = ' Nonlocal part applied using real spherical harmonics '
760 : : case default
761 [ + + - ]: 2287 : ABI_ERROR(sjoin("Invalid useylm: ", itoa(psps%useylm)))
762 : : end select
763 : 2287 : call wrtout(units,msg,mode)
764 : :
765 : 2287 : write(msg,'(a,i3)')' Max number of non-local projectors over l and type ',psps%mproj
766 : 2287 : call wrtout(units,msg,mode)
767 : :
768 : : write(msg,'(a,i3,2a,i3,2a,i3)')&
769 : 2287 : ' Highest angular momentum +1 ....... ',psps%mpsang,ch10,&
770 : 2287 : ' Max number of (l,n) components .. ',psps%lnmax, ch10,&
771 : 4574 : ' Max number of (l,m,n) components .. ',psps%lmnmax
772 : 2287 : call wrtout(units,msg,mode)
773 : :
774 : : !FIXME for paw n1xccc==1
775 : : ! Non-linear Core correction
776 [ + + ]: 2287 : if (psps%n1xccc/=0) then
777 : 1388 : write(msg,'(3a,2(a,i4,a),2a)')ch10,&
778 : 1388 : ' Pseudo-Core Charge Info: ',ch10,&
779 : 1388 : ' Number of radial points for pseudo-core charge .. ',psps%n1xccc,ch10,&
780 : 1388 : ' XC core-correction treatment (optnlxccc) ........ ',psps%optnlxccc,ch10,&
781 : 2776 : ' Radius for pseudo-core charge for each type ..... ',ch10
782 : 1388 : call wrtout(units,msg,mode)
783 [ + + ]: 3267 : do itypat=1,psps%ntypat
784 : 1879 : write(msg,'(a,i4,a,f12.4)')' - Atom type ',itypat,' has pseudo-core radius .. ',psps%xcccrc(itypat)
785 : 3267 : call wrtout(units,msg,mode)
786 : : end do
787 : : end if
788 : :
789 : : ! Alchemical mixing
790 [ + + ]: 2287 : if (psps%mtypalch/=0) then
791 : 26 : write(msg,'(3a,3(a,i4,a))')ch10,&
792 : 26 : ' Calculation with alchemical mixing:',ch10,&
793 : 26 : ' Number of pure pseudoatoms .... ',psps%ntyppure,ch10,&
794 : 26 : ' Number of pseudos for mixing .. ',psps%npspalch,ch10,&
795 : 52 : ' Alchemical pseudoatoms ........ ',psps%ntypalch,ch10
796 : 26 : call wrtout(units,msg,mode)
797 [ + + ]: 78 : do ipsp_alch=1,psps%npspalch
798 [ + + ]: 130 : do ityp_alch=1,psps%ntypalch
799 : 104 : write(std_out,*)' mixalch ',psps%mixalch(ipsp_alch,ityp_alch)
800 : : end do
801 : : end do
802 [ + + ]: 52 : do ityp_alch=1,psps%ntypalch
803 : 26 : write(msg,'(a,i4,a,i4)')' For alchemical atom no. ',ityp_alch,' algalch is .. ',psps%algalch(ityp_alch)
804 : 52 : call wrtout(units,msg,mode)
805 : : end do
806 : : end if
807 : :
808 : : ! Info in Q-grid for spline of form factors
809 : 2287 : write(msg,'(3a,a,i6,a,a,i6)')ch10,&
810 : 2287 : ' Info on the Q-grid used for form factors in spline form: ',ch10,&
811 : 2287 : ' Number of q-points for radial functions ffspl .. ',psps%mqgrid_ff,ch10,&
812 : 4574 : ' Number of q-points for vlspl ................... ',psps%mqgrid_vl
813 : 2287 : call wrtout(units,msg,mode)
814 : :
815 [ + - ]: 2287 : if (psps%vlspl_recipSpace) then
816 : 2287 : call wrtout(units,' vloc is computed in Reciprocal Space ',mode)
817 : : else
818 : 0 : call wrtout(units,' vloc is computed in Real Space ',mode)
819 : : end if
820 [ + + ]: 2287 : if (psps%usepaw == 0) then
821 [ + + ]: 1746 : if (psps%nc_xccc_gspace == 0) call wrtout(units,' model core charge treated in real-space', mode)
822 [ + + ]: 1746 : if (psps%nc_xccc_gspace == 1) call wrtout(units,' model core charge treated in G-space', mode)
823 : : end if
824 : :
825 : : !TODO additional stuff that might be printed
826 : 2287 : call wrtout(units, "", mode)
827 [ + + ]: 5221 : do itypat=1,psps%ntypat
828 : 2934 : write(msg,'(a,i0,a,i0)')' XC functional for type ',itypat,' is ',psps%pspxc(itypat)
829 : 2934 : call wrtout(units,msg,mode)
830 : : !write(std_out,*)psps%ziontypat(itypat),psps%znucltypat(itypat)
831 [ + + ]: 5221 : if (psps%usepaw == 0) then
832 : 2163 : call wrtout(units, sjoin(" Pseudo valence available: ", yesno(psps%nctab(itypat)%has_tvale)), mode)
833 : : end if
834 : : end do
835 : :
836 : : !integer, pointer :: pspxc(:)
837 : : ! pspxc(ntypat)
838 : : ! For each type of psp, the XC functional that was used to generate it, as given by the psp file
839 [ + + ]: 2287 : if (my_prtvol>=3) then
840 [ + + ]: 463 : do ips=1,psps%npsp
841 [ + - ]: 256 : write(std_out,*)' Pseudo number ',ips,' read from ',trim(psps%filpsp(ips))
842 : 256 : write(std_out,*)' Format or code ',psps%pspcod(ips)
843 : 256 : write(std_out,*)' Generation date ',psps%pspdat(ips)
844 [ + - ]: 463 : write(std_out,*)' Content of first line: ', trim(psps%title(ips))
845 : : end do
846 : : end if
847 : :
848 : 2287 : call wrtout(units, "", mode)
849 : :
850 : 2287 : end subroutine psps_print
851 : : !!***
852 : :
853 : : !----------------------------------------------------------------------
854 : :
855 : : !!****f* m_psps/psps_ncwrite_path
856 : : !! NAME
857 : : !! psps_ncwrite_path
858 : : !!
859 : : !! FUNCTION
860 : : !! Create a new NETCDF file,
861 : : !! and output the most important arrays defined in the pseudopotential_type
862 : : !! for futher post-processing.
863 : : !! This function should be called by master node only.
864 : : !!
865 : : !! INPUTS
866 : : !! path=File name.
867 : : !!
868 : : !! SOURCE
869 : :
870 : 1 : subroutine psps_ncwrite_path(psps, path)
871 : :
872 : : !Arguments ------------------------------------
873 : : class(pseudopotential_type),intent(in) :: psps
874 : : character(len=*),intent(in) :: path
875 : :
876 : : !Local variables-------------------------------
877 : : integer :: ncid
878 : : ! *************************************************************************
879 : :
880 [ - + ]: 1 : NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self))
881 : 1 : call psps_ncwrite(psps, ncid)
882 [ - + ]: 1 : NCF_CHECK(nf90_close(ncid))
883 : :
884 : 1 : end subroutine psps_ncwrite_path
885 : : !!***
886 : :
887 : : !----------------------------------------------------------------------
888 : :
889 : : !!****f* m_psps/psps_ncwrite
890 : : !! NAME
891 : : !! psps_ncwrite
892 : : !!
893 : : !! FUNCTION
894 : : !! Output the most important arrays defined in the pseudopotential_type
895 : : !! in NETCDF file format for futher post-processing.
896 : : !! This function should be called by master node only.
897 : : !!
898 : : !! INPUTS
899 : : !! ncid=NC file handle.
900 : : !!
901 : : !! SOURCE
902 : :
903 : 3312 : subroutine psps_ncwrite(psps, ncid)
904 : :
905 : : !Arguments ------------------------------------
906 : : class(pseudopotential_type),intent(in) :: psps
907 : : integer,intent(in) :: ncid
908 : :
909 : : !Local variables-------------------------------
910 : : !scalars
911 : : integer :: ipsp,itypat,ncerr
912 : : integer :: with_xccc, n1xccc, with_alch
913 : : integer :: with_xcctau
914 : : !arrays
915 : 3312 : real(dp), allocatable :: dummy3(:,:,:)
916 : : !real(dp), allocatable :: dummy1(:)
917 : : ! *************************************************************************
918 : :
919 : 3312 : with_alch = 0 ! Alchemical IO not supported at the moment.
920 : : !psps%mtypalch = zero
921 : :
922 : : ! GA: Note that lnmax is not used in the DDB text format,
923 : : ! so lnmax and lmnmax may be inconsistent in the netcdf file.
924 : : !NCF_CHECK(nctk_set_defmode(ncid))
925 : :
926 : 3312 : with_xccc = 0
927 [ + + ]: 3312 : if (psps%n1xccc > 0) then
928 : 1943 : with_xccc = 1
929 : 1943 : with_xcctau = 1
930 : : end if
931 : 3312 : n1xccc = max(1, psps%n1xccc)
932 : :
933 [ + + + - : 3312 : if (.not. allocated(psps%xcccrc) .or. .not. allocated(psps%xccc1d) .or. psps%usepaw /= 0) then
+ + ]
934 : 1839 : with_xccc = 0
935 : : end if
936 [ + + + + ]: 3312 : if (.not. allocated(psps%xcctau1d) .or. psps%usepaw /= 0) then
937 : 1839 : with_xcctau = 0
938 : : end if
939 : :
940 : : ! Define dimensions
941 : : ncerr = nctk_def_dims(ncid, [ &
942 : : nctkdim_t("fnlen", fnlen + 1), &
943 : : nctkdim_t("md5_slen", md5_slen + 1), &
944 : : nctkdim_t("ntypat", psps%ntypat), &
945 : : nctkdim_t("npsp", psps%npsp), &
946 : : nctkdim_t("lnmax", psps%lnmax), &
947 : : nctkdim_t("lmnmax", psps%lmnmax), &
948 : : nctkdim_t("dimekb", psps%dimekb), &
949 : : nctkdim_t("mqgrid_vl", psps%mqgrid_vl), &
950 : : nctkdim_t("mqgrid_ff", psps%mqgrid_ff), &
951 : : nctkdim_t("n1xccc", n1xccc) &
952 [ + + ]: 36432 : ])
953 [ - + ]: 3312 : NCF_CHECK(ncerr)
954 : :
955 : : ! Define variables
956 : : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
957 [ + + ]: 19872 : "usepaw", "useylm", "with_xccc", "with_xcctau", "with_alch"])
958 [ - + ]: 3312 : NCF_CHECK(ncerr)
959 : :
960 : : ! Arrays
961 : : ncerr = nctk_def_arrays(ncid, [&
962 : : nctkarr_t("ziontypat", "dp", "ntypat"), &
963 : : nctkarr_t("znucltypat", "dp", "ntypat"), &
964 : : nctkarr_t("spinorbit", "int", "npsp"), &
965 : : nctkarr_t("qgrid_vl", "dp", "mqgrid_vl"), &
966 : : nctkarr_t("qgrid_ff", "dp", "mqgrid_ff"), &
967 : : nctkarr_t("vlspl", "dp", "mqgrid_vl, two, ntypat"), &
968 : : nctkarr_t("indlmn", "int", "six, lmnmax, ntypat"), &
969 : : nctkarr_t("ffspl", "dp", "mqgrid_ff, two, lnmax, ntypat"), &
970 : : nctkarr_t("filpsp", "char", "fnlen, npsp"), &
971 : : nctkarr_t("md5_pseudos", "char", "md5_slen, npsp") &
972 [ + + ]: 36432 : ])
973 [ - + ]: 3312 : NCF_CHECK(ncerr)
974 : :
975 [ + + ]: 3312 : if (psps%usepaw == 0) then
976 [ - + ]: 2685 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("ekb", "dp", "dimekb, ntypat")))
977 : : !if (with_xccc > 0) then
978 [ - + ]: 2685 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("xccc1d", "dp", "n1xccc, six, ntypat")))
979 [ + + ]: 2685 : if (with_xcctau > 0) then
980 [ - + ]: 1473 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("xcctau1d", "dp", "n1xccc, six, ntypat")))
981 : : end if
982 [ - + ]: 2685 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("xcccrc", "dp", "ntypat")))
983 : :
984 : : ncerr = nctk_def_arrays(ncid, [&
985 : : nctkarr_t("nc_tvalespl", "dp", "mqgrid_vl, two, ntypat"), &
986 : : nctkarr_t("nc_tcorespl", "dp", "mqgrid_vl, two, ntypat"), &
987 : : nctkarr_t("nc_ttaucorespl", "dp", "mqgrid_vl, two, ntypat") &
988 [ + + ]: 10740 : ])
989 [ - + ]: 2685 : NCF_CHECK(ncerr)
990 : : end if
991 : :
992 : : ! Write data
993 [ - + ]: 3312 : NCF_CHECK(nf90_put_var(ncid, vid("ziontypat"), psps%ziontypat))
994 [ - + ]: 3312 : NCF_CHECK(nf90_put_var(ncid, vid("znucltypat"), psps%znucltypat))
995 : : ! Note that znuclpsp and ziopsp are not read, since we set with_alch=0
996 : :
997 : : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
998 : : "usepaw", "useylm", "with_xccc", "with_xcctau", "with_alch"], &
999 [ + + + + ]: 36432 : [psps%usepaw, psps%useylm, with_xccc, with_xcctau, with_alch])
1000 [ - + ]: 3312 : NCF_CHECK(ncerr)
1001 : :
1002 [ + - ]: 3312 : if (allocated(psps%pspso)) then
1003 [ - + ]: 3312 : NCF_CHECK(nf90_put_var(ncid, vid("spinorbit"), psps%pspso))
1004 : : end if
1005 [ + + ]: 7531 : do ipsp=1,psps%npsp
1006 [ + + + + : 12657 : NCF_CHECK(nf90_put_var(ncid, vid("filpsp"), trim(psps%filpsp(ipsp)), start=[1, ipsp]))
- + - - -
- ]
1007 [ + + + + : 15969 : NCF_CHECK(nf90_put_var(ncid, vid("md5_pseudos"), trim(psps%md5_pseudos(ipsp)), start=[1, ipsp]))
- + - - -
- ]
1008 : : end do
1009 [ + + ]: 3312 : if (allocated(psps%qgrid_vl)) then
1010 [ - + ]: 3310 : NCF_CHECK(nf90_put_var(ncid, vid("qgrid_vl"), psps%qgrid_vl))
1011 : : end if
1012 [ + + ]: 3312 : if (allocated(psps%qgrid_ff)) then
1013 [ - + ]: 3310 : NCF_CHECK(nf90_put_var(ncid, vid("qgrid_ff"), psps%qgrid_ff))
1014 : : end if
1015 [ + - ]: 3312 : if (allocated(psps%indlmn)) then
1016 [ - + ]: 3312 : NCF_CHECK(nf90_put_var(ncid, vid("indlmn"), psps%indlmn))
1017 : : end if
1018 : :
1019 : : ! Local part in q-space and second derivative
1020 [ + + ]: 3312 : if (allocated(psps%vlspl)) then
1021 [ - + ]: 3310 : NCF_CHECK(nf90_put_var(ncid, vid("vlspl"), psps%vlspl))
1022 : : end if
1023 : :
1024 : : ! Form factors for each type of atom
1025 : : ! for each type and each (l,n) channel, ffnl(q) and second derivative
1026 [ + + ]: 3312 : if (allocated(psps%ffspl)) then
1027 [ - + ]: 3310 : NCF_CHECK(nf90_put_var(ncid, vid("ffspl"), psps%ffspl))
1028 : : end if
1029 : :
1030 [ + + ]: 3312 : if (with_xccc > 0) then
1031 : :
1032 : : ! Pseudo-core charge for each type of atom, on the real-space radial
1033 [ - + ]: 1473 : NCF_CHECK(nf90_put_var(ncid, vid("xcccrc"), psps%xcccrc))
1034 [ - + ]: 1473 : NCF_CHECK(nf90_put_var(ncid, vid("xccc1d"), psps%xccc1d))
1035 [ + - ]: 1473 : if (with_xcctau > 0) then
1036 [ - + ]: 1473 : NCF_CHECK(nf90_put_var(ncid, vid("xcctau1d"), psps%xcctau1d))
1037 : : end if
1038 : :
1039 : : !else
1040 : :
1041 : : ! ABI_MALLOC(dummy1, (psps%ntypat))
1042 : : ! dummy1 = zero
1043 : : ! NCF_CHECK(nf90_put_var(ncid, vid("xcccrc"), dummy1))
1044 : : ! ABI_FREE(dummy1)
1045 : :
1046 : : ! ABI_MALLOC(dummy3, (n1xccc, 6, psps%ntypat))
1047 : : ! dummy3 = zero
1048 : : ! NCF_CHECK(nf90_put_var(ncid, vid("xccc1d"), dummy3))
1049 : : ! ABI_FREE(dummy3)
1050 : :
1051 : : end if
1052 : :
1053 : : ! NC-only: add tcore_spl and tvalespl in q-space
1054 [ + + - + ]: 5997 : if (psps%usepaw == 0) then
1055 [ + - ]: 2685 : if (allocated(psps%ekb)) then
1056 [ - + ]: 2685 : NCF_CHECK(nf90_put_var(ncid, vid("ekb"), psps%ekb))
1057 : : end if
1058 [ + + ]: 5907 : do itypat=1,psps%ntypat
1059 : :
1060 : : ! TODO Could write variables has_tvale and has_tcore
1061 [ + + ]: 3222 : if (psps%nctab(itypat)%has_tvale) then
1062 [ + + ]: 2296 : ncerr = nf90_put_var(ncid, vid("nc_tvalespl"), psps%nctab(itypat)%tvalespl, start=[1,1,itypat])
1063 [ - + ]: 574 : NCF_CHECK(ncerr)
1064 : : else
1065 [ + + + - : 10592 : ABI_MALLOC(dummy3, (psps%mqgrid_vl, 2, psps%ntypat))
- + + - +
- + + - +
- + - + ]
1066 [ + + + + : 21551370 : dummy3 = zero
+ + ]
1067 : 5296 : ncerr = nf90_put_var(ncid, vid("nc_tvalespl"), dummy3)
1068 [ - + ]: 2648 : NCF_CHECK(ncerr)
1069 [ - + ]: 2648 : ABI_FREE(dummy3)
1070 : : end if
1071 [ + + ]: 5907 : if (psps%nctab(itypat)%has_tcore) then
1072 [ + + ]: 6880 : ncerr = nf90_put_var(ncid, vid("nc_tcorespl"), psps%nctab(itypat)%tcorespl, start=[1,1,itypat])
1073 [ - + ]: 1720 : NCF_CHECK(ncerr)
1074 [ + - ]: 1720 : if (with_xcctau > 0) then
1075 [ + + ]: 6880 : ncerr = nf90_put_var(ncid, vid("nc_ttaucorespl"), psps%nctab(itypat)%ttaucorespl, start=[1,1,itypat])
1076 [ - + ]: 1720 : NCF_CHECK(ncerr)
1077 : : end if
1078 : : else
1079 [ + + + - : 6008 : ABI_MALLOC(dummy3, (psps%mqgrid_vl, 2, psps%ntypat))
- + + - +
- + + - +
- + - + ]
1080 [ + + + + : 12713392 : dummy3 = zero
+ + ]
1081 : 3004 : ncerr = nf90_put_var(ncid, vid("nc_tcorespl"), dummy3)
1082 [ - + ]: 1502 : NCF_CHECK(ncerr)
1083 [ + + ]: 1502 : if (with_xcctau > 0) then
1084 : 182 : ncerr = nf90_put_var(ncid, vid("nc_ttaucorespl"), dummy3)
1085 [ - + ]: 91 : NCF_CHECK(ncerr)
1086 : : end if
1087 [ - + ]: 1502 : ABI_FREE(dummy3)
1088 : : end if
1089 : : end do
1090 : : end if
1091 : :
1092 : : contains
1093 : 50285 : integer function vid(vname)
1094 : : character(len=*),intent(in) :: vname
1095 : 50285 : vid = nctk_idname(ncid, vname)
1096 : : end function vid
1097 : :
1098 : : end subroutine psps_ncwrite
1099 : : !!***
1100 : :
1101 : : !----------------------------------------------------------------------
1102 : :
1103 : : !!****f* m_psps/psps_ncread
1104 : : !! NAME
1105 : : !! psps_ncread
1106 : : !!
1107 : : !! FUNCTION
1108 : : !! Read the most important arrays defined in the pseudopotential_type
1109 : : !! in NETCDF file format.
1110 : : !! This function should be called by master node only.
1111 : : !!
1112 : : !! INPUTS
1113 : : !!
1114 : : !! SOURCE
1115 : :
1116 : 563 : subroutine psps_ncread(psps, ncid)
1117 : :
1118 : : !Arguments ------------------------------------
1119 : : class(pseudopotential_type),intent(inout) :: psps
1120 : : integer,intent(in) :: ncid
1121 : :
1122 : : !Local variables-------------------------------
1123 : : !scalars
1124 : : integer :: ipsp,itypat, ncerr, with_xccc, with_xcctau
1125 : : ! *********************************************************************
1126 : :
1127 : : ! Note: Some dimensions and variables are written conditionally,
1128 : : ! so try to read those but ignore errors
1129 : 563 : call psps_free(psps)
1130 : :
1131 : 563 : psps%dimekb = zero
1132 : 563 : psps%lmnmax = zero
1133 : 563 : psps%lnmax = zero
1134 : 563 : psps%mproj = zero
1135 : 563 : psps%mpsang = zero
1136 : 563 : psps%mpspso = zero
1137 : 563 : psps%mpssoang = zero
1138 : 563 : psps%mqgrid_ff = zero
1139 : 563 : psps%mqgrid_vl = zero
1140 : 563 : psps%mtypalch = zero
1141 : 563 : psps%npsp = zero
1142 : 563 : psps%npspalch = zero
1143 : 563 : psps%ntypat = zero
1144 : 563 : psps%ntypalch = zero
1145 : 563 : psps%ntyppure = zero
1146 : 563 : psps%n1xccc = zero
1147 : 563 : psps%optnlxccc = zero
1148 : 563 : psps%positron = zero
1149 : 563 : psps%usepaw = zero
1150 : 563 : psps%usewvl = zero
1151 : 563 : psps%useylm = zero
1152 : 563 : psps%nc_xccc_gspace = zero
1153 : 563 : psps%vlspl_recipSpace = .false.
1154 : :
1155 : : ! Read dimensions
1156 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "ntypat", psps%ntypat))
1157 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "npsp", psps%npsp))
1158 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "lnmax", psps%lnmax))
1159 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "lmnmax", psps%lmnmax))
1160 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "dimekb", psps%dimekb))
1161 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "mqgrid_vl", psps%mqgrid_vl))
1162 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "mqgrid_ff", psps%mqgrid_ff))
1163 [ - + ]: 563 : NCF_CHECK(nctk_get_dim(ncid, "n1xccc", psps%n1xccc))
1164 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "usepaw"), psps%usepaw))
1165 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "useylm"), psps%useylm))
1166 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "with_xccc"), with_xccc))
1167 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "with_xcctau"), with_xcctau))
1168 : :
1169 [ + + ]: 563 : if (psps%usepaw > 0) then
1170 : 7 : with_xccc = 0
1171 : 7 : with_xcctau = 0
1172 : : end if
1173 [ + + ]: 563 : if (with_xccc == 0) psps%n1xccc = 0
1174 : :
1175 : : ! Allocate arrays
1176 : 563 : call psp2params_init(psps%gth_params, psps%npsp)
1177 [ - + + - : 1689 : ABI_MALLOC(psps%filpsp,(psps%npsp))
+ - + - -
+ - + -
+ ]
1178 [ - + + - : 1126 : ABI_MALLOC(psps%title,(psps%npsp))
- + - + -
+ ]
1179 [ - + + - : 1689 : ABI_MALLOC(psps%md5_pseudos,(psps%npsp))
+ - + - -
+ - + -
+ ]
1180 : :
1181 [ - + + - : 1689 : ABI_MALLOC(psps%pspcod,(psps%npsp))
+ - + - -
+ - + -
+ ]
1182 [ - + + - : 1126 : ABI_MALLOC(psps%pspdat,(psps%npsp))
- + - + -
+ ]
1183 [ - + + - : 1126 : ABI_MALLOC(psps%pspxc,(psps%npsp))
- + - + -
+ ]
1184 [ - + + - : 1126 : ABI_MALLOC(psps%pspso,(psps%npsp))
- + - + -
+ ]
1185 : :
1186 [ + + ]: 1331 : psps%pspcod = zero
1187 [ + + ]: 1331 : psps%pspdat = zero
1188 [ + + ]: 1331 : psps%pspxc = zero
1189 [ + + ]: 1331 : psps%pspso = zero
1190 : :
1191 : : ! GA: zionpsp and znuclpsp dont get written. We assume they are the same
1192 : : ! as ziontypat and znucltypat
1193 [ - + + - : 1689 : ABI_MALLOC(psps%zionpsp,(psps%npsp))
+ - + - -
+ - + -
+ ]
1194 [ - + + - : 1126 : ABI_MALLOC(psps%znuclpsp,(psps%npsp))
- + - + -
+ ]
1195 [ - + + - : 1689 : ABI_MALLOC(psps%ziontypat,(psps%ntypat))
+ - + - -
+ - + -
+ ]
1196 [ - + + - : 1126 : ABI_MALLOC(psps%znucltypat,(psps%ntypat))
- + - + -
+ ]
1197 [ - + + - : 1126 : ABI_MALLOC(psps%xcccrc,(psps%ntypat))
- + - + -
+ ]
1198 [ + + + - : 1689 : ABI_MALLOC(psps%qgrid_vl,(psps%mqgrid_vl))
+ - + + -
+ - + -
+ ]
1199 [ + + + - : 1689 : ABI_MALLOC(psps%qgrid_ff,(psps%mqgrid_ff))
+ - + + -
+ - + -
+ ]
1200 [ - + + - : 2252 : ABI_MALLOC(psps%indlmn,(6,psps%lmnmax,psps%ntypat))
- + + - +
- + - - +
- + - + ]
1201 [ + + + - : 2252 : ABI_MALLOC(psps%vlspl,(psps%mqgrid_vl,2,psps%ntypat))
- + + - +
- + + - +
- + - + ]
1202 [ + + + - : 2815 : ABI_MALLOC(psps%ffspl,(psps%mqgrid_ff,2,psps%lmnmax,psps%ntypat))
- + + - -
+ + - + -
+ + - + -
+ - + ]
1203 [ - + + - : 2252 : ABI_MALLOC(psps%ekb,(psps%dimekb,psps%ntypat * (1 - psps%usepaw)))
+ + + - +
- + + - +
- + - + ]
1204 [ + + + - : 2252 : ABI_MALLOC(psps%xccc1d,(psps%n1xccc,6,psps%ntypat))
- + + - +
- + + - +
- + - + ]
1205 [ + + + - : 1689 : ABI_MALLOC(psps%xcctau1d,(psps%n1xccc,6,psps%ntypat))
- + + - -
+ - + -
+ ]
1206 [ - + + - : 2457 : ABI_MALLOC(psps%nctab,(psps%ntypat))
+ - + - -
+ - + - +
+ + + - -
+ - + - +
- + - + -
+ - + - +
- + ]
1207 [ + + ]: 563 : if (psps%usepaw == 0) then
1208 [ + + ]: 1310 : do itypat=1,psps%ntypat
1209 : 754 : psps%nctab(itypat)%mqgrid_vl = psps%mqgrid_vl
1210 : 754 : psps%nctab(itypat)%dncdq0 = zero
1211 : 754 : psps%nctab(itypat)%d2ncdq0 = zero
1212 : 754 : psps%nctab(itypat)%dtaucdq0 = zero
1213 : 754 : psps%nctab(itypat)%d2taucdq0 = zero
1214 : 754 : psps%nctab(itypat)%dnvdq0 = zero
1215 : 754 : psps%nctab(itypat)%num_tphi = zero
1216 : 754 : psps%nctab(itypat)%has_jtot = .False.
1217 : :
1218 : 754 : psps%nctab(itypat)%has_tvale = .False.
1219 : 754 : psps%nctab(itypat)%has_tcore = .False.
1220 : : ! GA: Do we even need those?
1221 [ + + + - : 2262 : ABI_MALLOC(psps%nctab(itypat)%tvalespl,(psps%mqgrid_vl,2))
+ - + + -
+ - + -
+ ]
1222 [ + + + - : 2262 : ABI_MALLOC(psps%nctab(itypat)%tcorespl,(psps%mqgrid_vl,2))
+ - + + -
+ - + -
+ ]
1223 [ + + + - : 2262 : ABI_MALLOC(psps%nctab(itypat)%ttaucorespl,(psps%mqgrid_vl,2))
+ - + + -
+ - + -
+ ]
1224 [ + + + + ]: 4455898 : psps%nctab(itypat)%tvalespl = zero
1225 [ + + + + ]: 4455898 : psps%nctab(itypat)%tcorespl = zero
1226 [ + + + + ]: 4456454 : psps%nctab(itypat)%ttaucorespl = zero
1227 : : end do
1228 : : end if
1229 : :
1230 [ + + + + ]: 3990 : psps%ekb = zero
1231 [ + + + + : 25320 : psps%indlmn = zero
+ + ]
1232 [ + + + + : 7568963 : psps%xccc1d = zero
+ + ]
1233 [ + + + + : 7568963 : psps%xcctau1d = zero
+ + ]
1234 [ + + ]: 1331 : psps%xcccrc = zero
1235 [ + + + + : 4540531 : psps%vlspl = zero
+ + ]
1236 [ + + + + : 20004426 : psps%ffspl = zero
+ + + + ]
1237 [ + + ]: 1678198 : psps%qgrid_vl = zero
1238 [ + + ]: 1678198 : psps%qgrid_ff = zero
1239 : :
1240 [ + + ]: 1331 : psps%zionpsp = zero
1241 [ + + ]: 1331 : psps%znuclpsp = zero
1242 [ + + ]: 1331 : psps%ziontypat = zero
1243 [ + + ]: 1331 : psps%znucltypat = zero
1244 : :
1245 : : ! Read variables
1246 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "ziontypat"), psps%ziontypat))
1247 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "znucltypat"), psps%znucltypat))
1248 : : ! Not dealing with alchemical at the moment.
1249 : :
1250 : :
1251 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "spinorbit"), psps%pspso))
1252 : :
1253 [ + + ]: 1331 : do ipsp=1,psps%npsp
1254 [ + + ]: 2304 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "filpsp"), psps%filpsp(ipsp), start=[1,ipsp])
1255 [ + + ]: 2304 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "md5_pseudos"), psps%md5_pseudos(ipsp), start=[1,ipsp])
1256 : 1331 : psps%title(ipsp) = ''
1257 : : end do
1258 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "qgrid_vl"), psps%qgrid_vl))
1259 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "qgrid_ff"), psps%qgrid_ff))
1260 [ - + ]: 563 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "indlmn"), psps%indlmn))
1261 : 563 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "vlspl"), psps%vlspl)
1262 : 563 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "ffspl"), psps%ffspl)
1263 : :
1264 [ + + ]: 563 : if (psps%usepaw == 0) then
1265 : 556 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "ekb"), psps%ekb)
1266 : :
1267 [ + + ]: 556 : if (with_xccc > 0) then
1268 [ - + ]: 351 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "xcccrc"), psps%xcccrc))
1269 [ - + ]: 351 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "xccc1d"), psps%xccc1d))
1270 [ + - ]: 351 : if (with_xcctau > 0) then
1271 [ - + ]: 351 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "xcctau1d"), psps%xcctau1d))
1272 : : end if
1273 : : end if
1274 : :
1275 : : ! GA: Why bother reading it?
1276 [ + + ]: 1310 : do itypat=1,psps%ntypat
1277 [ + + ]: 3016 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "nc_tvalespl"), psps%nctab(itypat)%tvalespl, start=[1,1,itypat])
1278 [ + + ]: 3572 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "nc_tcorespl"), psps%nctab(itypat)%tcorespl, start=[1,1,itypat])
1279 : : end do
1280 [ + + ]: 556 : if (with_xcctau > 0) then
1281 [ + + ]: 855 : do itypat=1,psps%ntypat
1282 [ + + ]: 2367 : ncerr = nf90_get_var(ncid, nctk_idname(ncid, "nc_ttaucorespl"), psps%nctab(itypat)%ttaucorespl, start=[1,1,itypat])
1283 : : end do
1284 : : end if
1285 : :
1286 : : end if
1287 : :
1288 : 563 : end subroutine psps_ncread
1289 : : !!***
1290 : :
1291 : : !----------------------------------------------------------------------
1292 : :
1293 : : !!****f* m_psps/psp2params_init
1294 : : !! NAME
1295 : : !! psp2params_init
1296 : : !!
1297 : : !! FUNCTION
1298 : : !! Allocate and initialise the data structure holding parameters for the GTH
1299 : : !! pseudo-potentials.
1300 : : !!
1301 : : !! MJV note: this should be renamed: psp2 suggests it relates to pspcod 2,
1302 : : !! whereas it is actually 3
1303 : : !! the parameters would also be better off separated into C and h arrays
1304 : : !!
1305 : : !! INPUTS
1306 : : !! npsp=number of true pseudo used (not alchemy).
1307 : : !!
1308 : : !! OUTPUT
1309 : : !! gth_params <type (pseudopotential_gth_type)>=the values to allocate and initialise.
1310 : : !!
1311 : : !! SOURCE
1312 : :
1313 [ + - ]: 1973 : subroutine psp2params_init(gth_params, npsp)
1314 : :
1315 : : !Arguments ------------------------------------
1316 : : class(pseudopotential_gth_type),intent(out) :: gth_params
1317 : : integer,intent(in) :: npsp
1318 : : ! *********************************************************************
1319 : :
1320 : : !Check array, no params are currently set.
1321 [ - + + - : 5919 : ABI_MALLOC(gth_params%set,(npsp))
+ - + - -
+ - + -
+ ]
1322 [ + + ]: 4667 : gth_params%set(:) = .false.
1323 : :
1324 : : !Check array, have geometric information been filled?
1325 [ - + + - : 3946 : ABI_MALLOC(gth_params%hasGeometry,(npsp))
- + - + -
+ ]
1326 [ + + ]: 4667 : gth_params%hasGeometry(:) = .false.
1327 : :
1328 : : !Coefficients for local part and projectors
1329 [ - + + - : 5919 : ABI_MALLOC(gth_params%psppar,(0:4, 0:6, npsp))
+ - + - -
+ - + -
+ ]
1330 [ + + + + : 117815 : gth_params%psppar = zero
+ + ]
1331 : :
1332 : : !Coefficients for spin orbit part
1333 [ - + + - : 5919 : ABI_MALLOC(gth_params%psp_k_par,(1:4, 1:3, npsp))
+ - + - -
+ - + -
+ ]
1334 [ + + + + : 45077 : gth_params%psp_k_par = zero
+ + ]
1335 : :
1336 : : !Different radii
1337 [ - + + - : 5919 : ABI_MALLOC(gth_params%radii_cf,(npsp, 3))
+ - + - -
+ - + -
+ ]
1338 [ + + + + ]: 15974 : gth_params%radii_cf = zero
1339 : :
1340 : 1973 : end subroutine psp2params_init
1341 : : !!***
1342 : :
1343 : : !----------------------------------------------------------------------
1344 : :
1345 : : !!****f* m_psps/psp2params_copy
1346 : : !! NAME
1347 : : !! psp2params_copy
1348 : : !!
1349 : : !! FUNCTION
1350 : : !!
1351 : : !! INPUTS
1352 : : !!
1353 : : !! OUTPUT
1354 : : !!
1355 : : !! SOURCE
1356 : :
1357 : 3477 : subroutine psp2params_copy(gth_paramsin, gth_paramsout)
1358 : :
1359 : : !Arguments ------------------------------------
1360 : : class(pseudopotential_gth_type),intent(in) :: gth_paramsin
1361 : : class(pseudopotential_gth_type),intent(inout) :: gth_paramsout
1362 : : ! *********************************************************************
1363 : :
1364 [ + + ]: 3477 : if (allocated(gth_paramsin%psppar)) then
1365 [ - + ]: 3467 : call alloc_copy( gth_paramsin%psppar, gth_paramsout%psppar)
1366 : : end if
1367 [ + + ]: 3477 : if (allocated(gth_paramsin%radii_cf)) then
1368 [ - + ]: 3467 : call alloc_copy( gth_paramsin%radii_cf, gth_paramsout%radii_cf)
1369 : : end if
1370 [ + + ]: 3477 : if (allocated(gth_paramsin%psp_k_par)) then
1371 [ - + ]: 3467 : call alloc_copy( gth_paramsin%psp_k_par, gth_paramsout%psp_k_par)
1372 : : end if
1373 [ + + ]: 3477 : if (allocated(gth_paramsin%hasGeometry)) then
1374 [ - + ]: 3467 : call alloc_copy( gth_paramsin%hasGeometry, gth_paramsout%hasGeometry)
1375 : : end if
1376 [ + + ]: 3477 : if (allocated(gth_paramsin%set)) then
1377 [ - + ]: 3467 : call alloc_copy( gth_paramsin%set, gth_paramsout%set)
1378 : : end if
1379 : :
1380 : 3477 : end subroutine psp2params_copy
1381 : : !!***
1382 : :
1383 : : !----------------------------------------------------------------------
1384 : :
1385 : : !!****f* m_psps/psp2params_free
1386 : : !! NAME
1387 : : !! psp2params_free
1388 : : !!
1389 : : !! FUNCTION
1390 : : !! Deallocate a previously allocated data structure for storage of GTH parameters.
1391 : : !!
1392 : : !! INPUTS
1393 : : !!
1394 : : !! SIDE EFFECTS
1395 : : !! gth_params <type (pseudopotential_gth_type)>=the values to deallocate.
1396 : : !!
1397 : : !! SOURCE
1398 : :
1399 : 7051 : subroutine psp2params_free(gth_params)
1400 : :
1401 : : !Arguments ------------------------------------
1402 : : class(pseudopotential_gth_type),intent(inout) :: gth_params
1403 : : ! *********************************************************************
1404 : :
1405 [ + + ]: 7051 : ABI_SFREE(gth_params%set)
1406 [ + + ]: 7051 : ABI_SFREE(gth_params%hasGeometry)
1407 : :
1408 : : ! Coefficients for local part and projectors
1409 [ + + ]: 7051 : ABI_SFREE(gth_params%psppar)
1410 : :
1411 : : ! Coefficients for spin orbit part
1412 [ + + ]: 7051 : ABI_SFREE(gth_params%psp_k_par)
1413 : :
1414 : : ! Different radii
1415 [ + + ]: 7051 : ABI_SFREE(gth_params%radii_cf)
1416 : :
1417 : 7051 : end subroutine psp2params_free
1418 : : !!***
1419 : :
1420 : : !!****f* m_psps/nctab_init
1421 : : !! NAME
1422 : : !! nctab_init
1423 : : !!
1424 : : !! FUNCTION
1425 : : !! Create nctab_t.
1426 : : !!
1427 : : !! INPUTS
1428 : : !! mqgrid_vl=Number of q-points
1429 : : !! has_tcore=True if the pseudo has NLCC.
1430 : : !! has_tvale=True if the atomic valence density is available.
1431 : : !!
1432 : : !! SOURCE
1433 : :
1434 : 1411 : subroutine nctab_init(nctab, mqgrid_vl, has_tcore, has_tvale)
1435 : :
1436 : : !Arguments ------------------------------------
1437 : : class(nctab_t),intent(inout) :: nctab
1438 : : integer,intent(in) :: mqgrid_vl
1439 : : logical,intent(in) :: has_tcore, has_tvale
1440 : : ! *************************************************************************
1441 : :
1442 : 1411 : nctab%mqgrid_vl = mqgrid_vl
1443 : :
1444 : : ! The array for the model core charge is always allocated and initialized with zeros.
1445 : : ! This approach is similar to the one used in the PAW code.
1446 : : ! has_tcore tells us whether the model core charge is present or not.
1447 : 1411 : nctab%has_tcore = has_tcore
1448 : 1411 : nctab%dncdq0 = zero; nctab%d2ncdq0 = zero
1449 : 1411 : nctab%dtaucdq0 = zero; nctab%d2taucdq0 = zero
1450 [ - + + - : 8683775 : ABI_CALLOC(nctab%tcorespl, (mqgrid_vl, 2))
+ - + - -
+ - + - +
+ + + + ]
1451 [ - + + - : 8682364 : ABI_CALLOC(nctab%ttaucorespl, (mqgrid_vl, 2))
- + - + -
+ + + +
+ ]
1452 : :
1453 : : ! tvalespl is allocated only if available.
1454 : 1411 : nctab%has_tvale = has_tvale
1455 : 1411 : nctab%dnvdq0 = zero
1456 [ - + ]: 1411 : if (has_tvale) then
1457 [ # # # # : 0 : ABI_CALLOC(nctab%tvalespl, (mqgrid_vl, 2))
# # # # #
# # # #
# ]
1458 : : end if
1459 : :
1460 : 1411 : end subroutine nctab_init
1461 : : !!***
1462 : :
1463 : : !!****f* m_psps/nctab_free
1464 : : !! NAME
1465 : : !! nctab_free
1466 : : !!
1467 : : !! FUNCTION
1468 : : !! Free memory allocated in nctab_t
1469 : : !!
1470 : : !! SOURCE
1471 : :
1472 : 6291 : subroutine nctab_free(nctab)
1473 : :
1474 : : !Arguments ------------------------------------
1475 : : class(nctab_t),intent(inout) :: nctab
1476 : : ! *************************************************************************
1477 : :
1478 [ + + ]: 6291 : ABI_SFREE(nctab%tvalespl)
1479 [ + + ]: 6291 : ABI_SFREE(nctab%tvaletauspl)
1480 [ + + ]: 6291 : ABI_SFREE(nctab%tcorespl)
1481 [ + + ]: 6291 : ABI_SFREE(nctab%ttaucorespl)
1482 [ + + ]: 6291 : ABI_SFREE(nctab%tphi_qspl)
1483 [ + + ]: 6291 : ABI_SFREE(nctab%tphi_n)
1484 [ + + ]: 6291 : ABI_SFREE(nctab%tphi_l)
1485 [ + + ]: 6291 : ABI_SFREE(nctab%tphi_jtot)
1486 [ + + ]: 6291 : ABI_SFREE(nctab%tphi_occ)
1487 : :
1488 : 6291 : end subroutine nctab_free
1489 : : !!***
1490 : :
1491 : : !!****f* m_psps/nctab_copy
1492 : : !! NAME
1493 : : !! nctab_copy
1494 : : !!
1495 : : !! FUNCTION
1496 : : !! Copy the object.
1497 : : !!
1498 : : !! SOURCE
1499 : :
1500 : 3371 : subroutine nctab_copy(nctabin, nctabout)
1501 : :
1502 : : !Arguments ------------------------------------
1503 : : class(nctab_t),intent(in) :: nctabin
1504 : : class(nctab_t),intent(inout) :: nctabout
1505 : : ! *************************************************************************
1506 : :
1507 : 3371 : nctabout%mqgrid_vl = nctabin%mqgrid_vl
1508 : 3371 : nctabout%has_tvale = nctabin%has_tvale
1509 : 3371 : nctabout%has_tcore = nctabin%has_tcore
1510 : 3371 : nctabout%dncdq0 = nctabin%dncdq0
1511 : 3371 : nctabout%d2ncdq0 = nctabin%d2ncdq0
1512 : 3371 : nctabout%dtaucdq0 = nctabin%dtaucdq0
1513 : 3371 : nctabout%d2taucdq0 = nctabin%d2taucdq0
1514 : 3371 : nctabout%dnvdq0 = nctabin%dnvdq0
1515 : 3371 : nctabout%has_tvaletau = nctabin%has_tvaletau
1516 : 3371 : nctabout%dnvtaudq0 = nctabin%dnvtaudq0
1517 : :
1518 : : ! TODO Why not check values of has_tvale and has_tcore?
1519 [ + + - + ]: 3371 : if (allocated(nctabin%tvalespl)) call alloc_copy(nctabin%tvalespl, nctabout%tvalespl)
1520 [ + + - + ]: 3371 : if (allocated(nctabin%tvaletauspl)) call alloc_copy(nctabin%tvaletauspl, nctabout%tvaletauspl)
1521 [ + + - + ]: 3371 : if (allocated(nctabin%tcorespl)) call alloc_copy(nctabin%tcorespl, nctabout%tcorespl)
1522 [ + + - + ]: 3371 : if (allocated(nctabin%ttaucorespl)) call alloc_copy(nctabin%ttaucorespl, nctabout%ttaucorespl)
1523 : :
1524 : 3371 : end subroutine nctab_copy
1525 : : !!***
1526 : :
1527 : : !!****f* m_psps/nctab_eval_tvalespl
1528 : : !! NAME
1529 : : !! nctab_eval_tvalespl
1530 : : !!
1531 : : !! FUNCTION
1532 : : !! Evalute spline-fit of the atomic pseudo valence charge in reciprocal space.
1533 : : !!
1534 : : !! INPUTS
1535 : : !! zion=nominal valence of atom as specified in psp file. Used to rescale the f(q=0) component
1536 : : !! mesh<pawrad_type>Radial mesh (r-space) used for the valence denity.
1537 : : !! valr(mesh%mesh_size)=Valence density in real space.
1538 : : !! mqgrid_vl=Number of points in the reciprocal space grid
1539 : : !! qgrid_vl(mqgrid_vl)=The coordinates of all the points of the radial q-grid
1540 : : !!
1541 : : !! SIDE EFFECTS
1542 : : !! nctabl%tvalspl(mqgrid_vl,2)
1543 : : !! nctab%dnvdq0
1544 : : !! nctab%d2nvdq0
1545 : : !!
1546 : : !! SOURCE
1547 : :
1548 : 503 : subroutine nctab_eval_tvalespl(nctab, zion, mesh, valr, mqgrid_vl, qgrid_vl)
1549 : :
1550 : : !Arguments ------------------------------------
1551 : : class(nctab_t),intent(inout) :: nctab
1552 : : integer,intent(in) :: mqgrid_vl
1553 : : real(dp),intent(in) :: zion
1554 : : type(pawrad_type),intent(in) :: mesh
1555 : : !arrays
1556 : : real(dp),intent(in) :: valr(mesh%mesh_size),qgrid_vl(mqgrid_vl)
1557 : :
1558 : : !Local variables-------------------------------
1559 : : real(dp) :: fact,yp1,ypn,d2nvdq0
1560 : : ! *************************************************************************
1561 : :
1562 : 503 : nctab%has_tvale = .True.
1563 [ + + ]: 503 : if (.not. allocated(nctab%tvalespl)) then
1564 [ - + + - : 1047 : ABI_MALLOC(nctab%tvalespl, (mqgrid_vl, 2))
+ - + - -
+ - + -
+ ]
1565 : : else
1566 [ - + ]: 154 : ABI_CHECK(size(nctab%tvalespl, dim=1) == mqgrid_vl, "wrong mqgrid_vl")
1567 : : end if
1568 : :
1569 : 503 : call pawpsp_cg(nctab%dnvdq0, d2nvdq0, mqgrid_vl, qgrid_vl, nctab%tvalespl(:,1), mesh, valr, yp1, ypn)
1570 [ + - + + ]: 295569 : call simp_gen(yp1, mesh%rad**2 * valr, mesh)
1571 : 503 : write(std_out,*)" valence charge (before rescaling) integrates to: ",four_pi*yp1
1572 : :
1573 : : ! Rescale the integral to have the correct number of valence electrons.
1574 : : ! In some cases, indeed, the radial mesh is not large enough and some valence charge is missing
1575 : : ! pawpsp_cg extrapolates the integrand beyond rmax but this is not enough.
1576 : : ! Remember that tvalespl is used to build an initial guess for rhor hence it's very important
1577 : : ! to have the correct electrostatic.
1578 : 503 : fact = zion / nctab%tvalespl(1,1)
1579 [ + + ]: 1528389 : nctab%tvalespl(:,1) = nctab%tvalespl(:,1) * fact
1580 : :
1581 : : ! Compute second derivative of tvalespl(q)
1582 : 503 : call paw_spline(qgrid_vl,nctab%tvalespl(:,1),mqgrid_vl,yp1,ypn,nctab%tvalespl(:,2))
1583 : :
1584 : 503 : end subroutine nctab_eval_tvalespl
1585 : : !!***
1586 : :
1587 : : !!****f* m_psps/nctab_eval_tvaletauspl
1588 : : !! NAME
1589 : : !! nctab_eval_tvaletauspl
1590 : : !!
1591 : : !! FUNCTION
1592 : : !! Evaluate spline-fit of the atomic pseudo valence kinetic energy density in reciprocal space.
1593 : : !!
1594 : : !! INPUTS
1595 : : !! mesh<pawrad_type>Radial mesh (r-space) used for the valence kinetic energy density.
1596 : : !! tauvalr(mesh%mesh_size)=Valence kinetic energy density in real space.
1597 : : !! mqgrid_vl=Number of points in the reciprocal space grid
1598 : : !! qgrid_vl(mqgrid_vl)=The coordinates of all the points of the radial q-grid
1599 : : !!
1600 : : !! SIDE EFFECTS
1601 : : !! nctab%tvaletauspl(mqgrid_vl,2)
1602 : : !! nctab%dnvtaudq0
1603 : : !!
1604 : : !! SOURCE
1605 : :
1606 : 1 : subroutine nctab_eval_tvaletauspl(nctab, mesh, tauvalr, mqgrid_vl, qgrid_vl)
1607 : :
1608 : : !Arguments ------------------------------------
1609 : : class(nctab_t),intent(inout) :: nctab
1610 : : integer,intent(in) :: mqgrid_vl
1611 : : type(pawrad_type),intent(in) :: mesh
1612 : : !arrays
1613 : : real(dp),intent(in) :: tauvalr(mesh%mesh_size),qgrid_vl(mqgrid_vl)
1614 : :
1615 : : !Local variables-------------------------------
1616 : : real(dp) :: yp1,ypn,d2nvtaudq0
1617 : : ! *************************************************************************
1618 : :
1619 : 1 : nctab%has_tvaletau = .True.
1620 [ + - ]: 1 : if (.not. allocated(nctab%tvaletauspl)) then
1621 [ - + + - : 3 : ABI_MALLOC(nctab%tvaletauspl, (mqgrid_vl, 2))
+ - + - -
+ - + -
+ ]
1622 : : else
1623 [ # # ]: 0 : ABI_CHECK(size(nctab%tvaletauspl, dim=1) == mqgrid_vl, "wrong mqgrid_vl")
1624 : : end if
1625 : :
1626 : 1 : call pawpsp_cg(nctab%dnvtaudq0, d2nvtaudq0, mqgrid_vl, qgrid_vl, nctab%tvaletauspl(:,1), mesh, tauvalr, yp1, ypn)
1627 : :
1628 : : ! No rescaling for kinetic energy density (unlike valence charge, we have no constraint on the integral).
1629 : :
1630 : : ! Compute second derivative of tvaletauspl(q)
1631 : 1 : call paw_spline(qgrid_vl,nctab%tvaletauspl(:,1),mqgrid_vl,yp1,ypn,nctab%tvaletauspl(:,2))
1632 : :
1633 : 1 : end subroutine nctab_eval_tvaletauspl
1634 : : !!***
1635 : :
1636 : : !!****f* m_psps/nctab_eval_tcorespl
1637 : : !! NAME
1638 : : !! nctab_eval_tcorespl
1639 : : !!
1640 : : !! FUNCTION
1641 : : !! Evalute spline-fit of the model core charge in reciprocal space.
1642 : : !!
1643 : : !! INPUTS
1644 : : !! xcccrc=maximum radius of the pseudo-core charge
1645 : : !! n1xccc=Number of radial points for the description of the pseudo-core charge
1646 : : !! (in the framework of the non-linear XC core correction)
1647 : : !! mqgrid_vl=Number of points in the reciprocal space grid
1648 : : !! qgrid_vl(mqgrid_vl)=The coordinates of all the points of the radial q-grid
1649 : : !! xccc1d(n1xccc,6)= The component xccc1d(n1xccc,1) is the pseudo-core charge
1650 : : !! on the radial grid. The components xccc1d(n1xccc,ideriv) give the ideriv-th derivative of the
1651 : : !! pseudo-core charge with respect to the radial distance.
1652 : : !! xcctau1d(n1xccc,6)= The component xcctau1d(n1xccc,1) is the pseudo-core kinetic energy density
1653 : : !! on the radial grid. The components xcctau1d(n1xccc,ideriv) give the ideriv-th derivative of the
1654 : : !! pseudo-core kinE den with respect to the radial distance.
1655 : : !!
1656 : : !! SIDE EFFECTS
1657 : : !! nctabl%tcorespl(mqgrid_vl,2)
1658 : : !! nctab%d2ncdq0
1659 : : !! nctab%dncdq0
1660 : : !!
1661 : : !! nctabl%ttaucorespl(mqgrid_vl,2)
1662 : : !! nctab%d2taucdq0
1663 : : !! nctab%dtaucdq0
1664 : : !!
1665 : : !! SOURCE
1666 : :
1667 : 2189 : subroutine nctab_eval_tcorespl(nctab, n1xccc, xcccrc, xccc1d, xcctau1d, mqgrid_vl, qgrid_vl)
1668 : :
1669 : : !Arguments ------------------------------------
1670 : : !scalars
1671 : : class(nctab_t),intent(inout) :: nctab
1672 : : integer,intent(in) :: n1xccc,mqgrid_vl
1673 : : real(dp),intent(in) :: xcccrc
1674 : : !arrays
1675 : : real(dp),intent(in) :: xccc1d(n1xccc,6),qgrid_vl(mqgrid_vl), xcctau1d(n1xccc,6)
1676 : :
1677 : : !Local variables-------------------------------
1678 : : real(dp) :: amesh,yp1,ypn
1679 : 2189 : type(pawrad_type) :: core_mesh
1680 : : ! *************************************************************************
1681 : :
1682 [ - + ]: 2189 : ABI_CHECK(mqgrid_vl == nctab%mqgrid_vl, "wrong mqgrid_vl")
1683 : :
1684 [ - + ]: 2189 : if (.not. allocated(nctab%tcorespl)) then
1685 [ # # # # : 0 : ABI_CALLOC(nctab%tcorespl, (mqgrid_vl, 2))
# # # # #
# # # # #
# # # # ]
1686 : : else
1687 [ - + ]: 2189 : ABI_CHECK(size(nctab%tcorespl, dim=1) == mqgrid_vl, "wrong mqgrid_vl")
1688 : : end if
1689 : :
1690 [ - + ]: 2189 : if (.not. allocated(nctab%ttaucorespl)) then
1691 [ # # # # : 0 : ABI_CALLOC(nctab%ttaucorespl, (mqgrid_vl, 2))
# # # # #
# # # # #
# # # # ]
1692 : : else
1693 [ - + ]: 2189 : ABI_CHECK(size(nctab%ttaucorespl, dim=1) == mqgrid_vl, "wrong mqgrid_vl")
1694 : : end if
1695 : :
1696 : : ! Skip loop if this atom has no core charge
1697 [ + + ]: 2189 : if (abs(xcccrc) < tol16) then
1698 : 1130 : nctab%has_tcore = .False.
1699 : : return
1700 : : end if
1701 : :
1702 : 1059 : nctab%has_tcore = .True.
1703 : : ! XCCC is given on a linear mesh.
1704 : 1059 : amesh = xcccrc / dble(n1xccc-1)
1705 : 1059 : call pawrad_init(core_mesh, mesh_size=n1xccc, mesh_type=1, rstep=amesh)
1706 : :
1707 : : ! Compute 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 n(r))dr].
1708 : : ! write(std_out,*)"xccc1d: amesh, min, max, minloc ",amesh,maxval(xccc1d(:,1)),minval(xccc1d(:,1)),minloc(xccc1d(:,1))
1709 : : call pawpsp_cg(nctab%dncdq0, nctab%d2ncdq0, mqgrid_vl, qgrid_vl, nctab%tcorespl(:,1), &
1710 : 1059 : core_mesh, xccc1d(:,1), yp1, ypn)
1711 : :
1712 : : ! Compute second derivative of tcorespl(q)
1713 : 1059 : call paw_spline(qgrid_vl, nctab%tcorespl(:,1), mqgrid_vl, yp1, ypn, nctab%tcorespl(:,2))
1714 : :
1715 : : ! idem for kinetic energy density
1716 : : call pawpsp_cg(nctab%dtaucdq0, nctab%d2taucdq0, mqgrid_vl, qgrid_vl, nctab%ttaucorespl(:,1), &
1717 : 1059 : core_mesh, xcctau1d(:,1), yp1, ypn)
1718 : 1059 : call paw_spline(qgrid_vl, nctab%ttaucorespl(:,1), mqgrid_vl, yp1, ypn, nctab%ttaucorespl(:,2))
1719 : :
1720 : 1059 : call pawrad_free(core_mesh)
1721 : :
1722 [ - + - + : 1059 : end subroutine nctab_eval_tcorespl
- + ]
1723 : : !!***
1724 : :
1725 : : !!****f* m_psps/nctab_mixalch
1726 : : !! NAME
1727 : : !! nctab_mixalch
1728 : : !!
1729 : : !! FUNCTION
1730 : : !! Mix the pseudopotential tables. Used for alchemical mixing.
1731 : : !!
1732 : : !! INPUTS
1733 : : !! nctabs(npspalch)=NC tables to be mixed
1734 : : !! npspalch=Number of alchemical pseudos.
1735 : : !! ntypalch=Number of types of alchemical pseudoatoms
1736 : : !! algalch(ntypalch)=For each type of pseudo atom, the algorithm to mix the pseudopotentials
1737 : : !! mixalch(npspalch,ntypalch)=Mixing coefficients to generate alchemical pseudo atoms
1738 : : !!
1739 : : !! OUTPUT
1740 : : !! mixtabs(ntypalch)=NC tables describing the alchemical pseudos
1741 : : !!
1742 : : !! SOURCE
1743 : :
1744 : 26 : subroutine nctab_mixalch(nctabs, npspalch, ntypalch, algalch, mixalch, mixtabs)
1745 : :
1746 : : !Arguments ------------------------------------
1747 : : !scalars
1748 : : integer,intent(in) :: npspalch,ntypalch
1749 : : !arrays
1750 : : integer,intent(in) :: algalch(ntypalch)
1751 : : real(dp),intent(in) :: mixalch(npspalch, ntypalch)
1752 : : type(nctab_t),intent(in) :: nctabs(npspalch)
1753 : : type(nctab_t),target,intent(inout) :: mixtabs(ntypalch)
1754 : :
1755 : : !Local variables-------------------------------
1756 : : !scalars
1757 : : integer :: ipspalch,itypalch
1758 : : logical :: has_tcore, has_tvale
1759 : : real(dp) :: mc
1760 : : type(nctab_t),pointer :: mix
1761 : : ! *************************************************************************
1762 : :
1763 [ + + + - : 78 : ABI_CHECK(all(nctabs(:)%mqgrid_vl == nctabs(1)%mqgrid_vl), "Wrong mqgrid_vl")
- + ]
1764 [ + + + - : 52 : ABI_CHECK(all(algalch == 1), "algalch /= 1 not implemented")
- + ]
1765 : :
1766 [ + + ]: 52 : do itypalch=1,ntypalch
1767 : :
1768 : : ! has_tcore is true is at least one pseudo has nlcc.
1769 : : ! has_tvale is true if *all* mixed pseudos have the PS valence charge.
1770 : 26 : has_tcore = .False.; has_tvale = .True.
1771 [ + + ]: 78 : do ipspalch=1,npspalch
1772 [ + + ]: 52 : if (abs(mixalch(ipspalch,itypalch)) < tol6) cycle
1773 [ + + ]: 42 : if (nctabs(ipspalch)%has_tcore) has_tcore = .True.
1774 [ + - ]: 68 : if (.not. nctabs(ipspalch)%has_tvale) has_tvale = .False.
1775 : : end do
1776 : : !write(std_out,*)has_tvale, has_tcore
1777 : :
1778 : 26 : call nctab_free(mixtabs(itypalch))
1779 : 26 : call nctab_init(mixtabs(itypalch), nctabs(1)%mqgrid_vl, has_tcore, has_tvale)
1780 : 26 : mix => mixtabs(itypalch)
1781 : :
1782 [ + + ]: 104 : do ipspalch=1,npspalch
1783 : 52 : mc = mixalch(ipspalch,itypalch)
1784 [ + + ]: 52 : if (abs(mc) < tol6) cycle
1785 : : ! Linear combination of the quantities
1786 : : ! Mix core for NLCC
1787 [ + + ]: 42 : if (has_tcore) then
1788 [ - + + + : 36036 : mix%tcorespl = mix%tcorespl + mc * nctabs(ipspalch)%tcorespl
+ + ]
1789 [ - + + + : 36036 : mix%ttaucorespl = mix%ttaucorespl + mc * nctabs(ipspalch)%ttaucorespl
+ + ]
1790 : 6 : mix%dncdq0 = mix%dncdq0 + mc * nctabs(ipspalch)%dncdq0
1791 : 6 : mix%d2ncdq0 = mix%d2ncdq0 + mc * nctabs(ipspalch)%d2ncdq0
1792 : 6 : mix%dtaucdq0 = mix%dtaucdq0 + mc * nctabs(ipspalch)%dtaucdq0
1793 : 6 : mix%d2taucdq0 = mix%d2taucdq0 + mc * nctabs(ipspalch)%d2taucdq0
1794 : : end if
1795 : : ! Mix pseudo valence charge.
1796 [ - + ]: 68 : if (has_tvale) then
1797 [ # # # # : 0 : mix%tvalespl = mix%tvalespl + mc * nctabs(ipspalch)%tvalespl
# # ]
1798 : 0 : mix%dnvdq0 = mix%dnvdq0 + mc * nctabs(ipspalch)%dnvdq0
1799 : : end if
1800 : : end do
1801 : :
1802 : : end do
1803 : :
1804 : 26 : end subroutine nctab_mixalch
1805 : : !!***
1806 : :
1807 : : end module m_psps
1808 : : !!***
|
