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 2718 : 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 2718 : usexml = 0; xmlpaw = 0; useupf = 0
113 :
114 2718 : if (open_file(path,msg,newunit=temp_unit,form='formatted',status='old') /= 0) then
115 0 : ABI_ERROR(msg)
116 : end if
117 2718 : rewind (unit=temp_unit,err=10,iomsg=errmsg)
118 :
119 2718 : read(temp_unit, "(a)",err=10,iomsg=errmsg) testxml
120 2718 : if(testxml(1:5)=='<?xml')then
121 321 : usexml = 1
122 321 : read(temp_unit,*,err=10,iomsg=errmsg) testxml
123 321 : if(testxml(1:4)=='<paw') xmlpaw = 1
124 : else
125 2397 : usexml = 0
126 2397 : 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 2718 : if (useupf == 0) then
144 2708 : rewind (unit=temp_unit,err=10,iomsg=errmsg)
145 2708 : read(temp_unit,*,err=10,iomsg=errmsg) testxml ! just a string, no relation to xml.
146 2708 : if(testxml(1:9)=='<PP_INFO>')then
147 1 : useupf = 1
148 : else
149 2707 : useupf = 0
150 : end if
151 : end if
152 :
153 2718 : close(unit=temp_unit,err=10,iomsg=errmsg)
154 :
155 2718 : 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 1392 : 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 4176 : ABI_MALLOC(psps%filpsp,(npsp))
196 4176 : ABI_MALLOC(psps%pspcod,(npsp))
197 2784 : ABI_MALLOC(psps%pspdat,(npsp))
198 2784 : ABI_MALLOC(psps%pspso,(npsp))
199 2784 : ABI_MALLOC(psps%pspxc,(npsp))
200 2784 : ABI_MALLOC(psps%title,(npsp))
201 4176 : ABI_MALLOC(psps%zionpsp,(npsp))
202 2784 : ABI_MALLOC(psps%znuclpsp,(npsp))
203 2784 : ABI_MALLOC(psps%epsatm,(npsp))
204 1392 : call psp2params_init(psps%gth_params, npsp)
205 :
206 3290 : psps%filpsp(1:npsp)=pspheads(1:npsp)%filpsp
207 3290 : psps%pspcod(1:npsp)=pspheads(1:npsp)%pspcod
208 3290 : psps%pspdat(1:npsp)=pspheads(1:npsp)%pspdat
209 3290 : psps%pspso(1:npsp)=pspheads(1:npsp)%pspso
210 3290 : psps%pspxc(1:npsp)=pspheads(1:npsp)%pspxc
211 3290 : psps%title(1:npsp)=pspheads(1:npsp)%title
212 3290 : psps%zionpsp(1:npsp)=pspheads(1:npsp)%zionpsp
213 3290 : psps%znuclpsp(1:npsp)=pspheads(1:npsp)%znuclpsp
214 :
215 : ! Transfer md5 checksum
216 4176 : ABI_MALLOC(psps%md5_pseudos, (npsp))
217 4682 : psps%md5_pseudos = pspheads(1:npsp)%md5_checksum
218 : !Set values independant from dtset
219 1392 : psps%npsp = npsp
220 : !Note that mpsang is the max of 1+lmax, with minimal value 1 (even for local psps, at present)
221 1392 : mpsang=1
222 1392 : n1xccc=pspheads(1)%xccc
223 3290 : do ii=1,psps%npsp
224 1898 : mpsang=max(pspheads(ii)%lmax+1,mpsang)
225 3290 : n1xccc=max(pspheads(ii)%xccc,n1xccc)
226 : end do
227 1392 : psps%mpsang = mpsang
228 1392 : psps%n1xccc = n1xccc
229 : ! Determine here whether the calculation is PAW
230 : ! If paw, all pspcod necessarily are 7 or 17 (see iofn2)
231 1392 : psps%usepaw =0
232 1392 : if (pspheads(1)%pspcod==7.or.pspheads(1)%pspcod==17) psps%usepaw=1
233 1392 : psps%mtypalch = mtypalch
234 :
235 1392 : 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 5285 : 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 5285 : 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 5285 : nptsgvec = 200 !This has to be chosen one and for all or else ??
280 5285 : newmqgrid = dtset%mqgrid
281 5285 : newmqgriddg = dtset%mqgriddg
282 :
283 : !JB:Which image to use ? I guess 1 always works
284 5285 : call matr3inv(dtset%rprimd_orig(:,:,1),gprimd_orig)
285 5285 : if ( dtset%usewvl == 0) then
286 : call setmqgrid(newmqgrid,newmqgriddg,dtset%ecut*dtset%dilatmx**2,&
287 5285 : 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 5285 : psps%mqgrid_ff = newmqgrid
292 5285 : if (psps%usepaw == 1) then
293 1176 : psps%mqgrid_vl = newmqgriddg
294 : else
295 4109 : 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 5285 : dtset%ntypat,dtset%ntypalch,pspheads)
301 :
302 5285 : psps%npspalch = dtset%npspalch
303 5285 : psps%ntypat = dtset%ntypat
304 5285 : psps%ntypalch = dtset%ntypalch
305 5285 : psps%ntyppure = dtset%ntyppure
306 :
307 : !Set the flag for reciprocal space or real space calculations
308 5285 : psps%vlspl_recipSpace = (dtset%icoulomb /= 1)
309 5285 : psps%positron = dtset%positron
310 5285 : psps%useylm = dtset%useylm
311 5285 : psps%usewvl = dtset%usewvl
312 :
313 : ! Define treatment of the model core density for NC pseudos.
314 5285 : psps%nc_xccc_gspace = dtset%nc_xccc_gspace
315 :
316 5285 : if (idtset > 1) then
317 3893 : ABI_SFREE(psps%algalch)
318 3893 : ABI_SFREE(psps%mixalch)
319 : end if
320 :
321 15855 : ABI_MALLOC(psps%algalch,(psps%ntypalch))
322 21140 : ABI_MALLOC(psps%mixalch,(psps%npspalch,psps%ntypalch))
323 5305 : psps%algalch(1:psps%ntypalch)=dtset%algalch(1:psps%ntypalch)
324 : !This value will be overwritten elsewhere in case there are different images ...
325 5345 : 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 5285 : psps%mpspso=1
331 12203 : do ipsp=1,dtset%npsp
332 12203 : if(dtset%nspinor==1)then
333 6403 : psps%pspso(ipsp)=0
334 :
335 : ! Ideally the following line should not exist, but at present, the space has to be booked
336 6403 : if(pspheads(ipsp)%pspso/=0)psps%mpspso=2
337 515 : else if (psps%usepaw==0) then
338 273 : if(dtset%so_psp(ipsp)/=1)then
339 126 : psps%pspso(ipsp)=dtset%so_psp(ipsp)
340 : else
341 147 : psps%pspso(ipsp)=pspheads(ipsp)%pspso
342 : end if
343 273 : if(psps%pspso(ipsp)/=0)psps%mpspso=2
344 273 : if(pspheads(ipsp)%pspso/=0)psps%mpspso=2
345 : else
346 242 : psps%pspso(ipsp)=1+dtset%pawspnorb
347 : end if
348 : end do
349 :
350 : !Set mpssoang, lmnmax, lnmax
351 5285 : if(psps%mpspso==1)then
352 4329 : psps%mpssoang=psps%mpsang
353 4329 : psps%lmnmax =lmnmax
354 4329 : psps%lnmax =lnmax
355 : else
356 956 : psps%mpssoang=2*psps%mpsang-1
357 956 : psps%lmnmax=lmnmaxso
358 956 : psps%lnmax=lnmaxso
359 : end if
360 :
361 : !T. Rangel: for wvl + paw do not change psps%lmnmax
362 5285 : if (psps%useylm==0 .and. psps%usepaw/=1 ) then
363 3973 : psps%lmnmax=psps%lnmax
364 : end if
365 :
366 : !Set dimekb
367 5285 : if (psps%usepaw==0) then
368 4109 : psps%dimekb=psps%lnmax
369 : else
370 1176 : 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 5285 : changed = 0
376 :
377 5285 : if(dimekb_old/=psps%dimekb .or. ntypat_old/=dtset%ntypat .or. usepaw_old/=psps%usepaw) then
378 1395 : changed = changed + 1
379 1395 : if(idtset/=1) then
380 3 : ABI_SFREE(psps%ekb)
381 : end if
382 5580 : ABI_MALLOC(psps%ekb,(psps%dimekb,dtset%ntypat*(1-psps%usepaw)))
383 8207 : psps%ekb = zero
384 1395 : dimekb_old=psps%dimekb
385 : end if
386 :
387 5285 : if(lmnmax_old/=psps%lmnmax .or. ntypat_old/=dtset%ntypat)then
388 1438 : changed = changed + 1
389 1438 : if(idtset/=1) then
390 46 : ABI_SFREE(psps%indlmn)
391 : end if
392 5752 : ABI_MALLOC(psps%indlmn,(6,psps%lmnmax,dtset%ntypat))
393 104024 : psps%indlmn = zero
394 1438 : lmnmax_old=psps%lmnmax
395 : end if
396 :
397 5285 : if(mqgridff_old/=psps%mqgrid_ff .or. lnmax_old/=psps%lnmax .or. ntypat_old/=dtset%ntypat)then
398 1397 : changed = changed + 1
399 1397 : if(idtset/=1) then
400 5 : ABI_SFREE(psps%ffspl)
401 5 : ABI_SFREE(psps%qgrid_ff)
402 : end if
403 6985 : ABI_MALLOC(psps%ffspl,(psps%mqgrid_ff,2,psps%lnmax,dtset%ntypat))
404 4191 : ABI_MALLOC(psps%qgrid_ff,(psps%mqgrid_ff))
405 49975914 : psps%ffspl = zero
406 4291039 : psps%qgrid_ff = zero
407 1397 : mqgridff_old=psps%mqgrid_ff
408 1397 : lnmax_old=psps%lnmax
409 : end if
410 :
411 5285 : if(mqgridvl_old/=psps%mqgrid_vl .or. ntypat_old/=dtset%ntypat)then
412 1395 : changed = changed + 1
413 1395 : 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 1395 : if (idtset/=1 .and. .not.psps%vlspl_recipSpace) then
424 0 : ABI_SFREE(psps%dvlspl)
425 : end if
426 :
427 4185 : ABI_MALLOC(psps%qgrid_vl,(psps%mqgrid_vl))
428 5580 : ABI_MALLOC(psps%vlspl,(psps%mqgrid_vl,2,dtset%ntypat))
429 4300663 : psps%qgrid_vl = zero
430 11644452 : psps%vlspl = zero
431 :
432 1395 : if (psps%usepaw == 0) then
433 : ! If you change usepaw in the input, you will get what you deserve!
434 4354 : ABI_MALLOC(psps%nctab, (dtset%ntypat))
435 2336 : do itypat=1,dtset%ntypat
436 2336 : call nctab_init(psps%nctab(itypat), psps%mqgrid_vl, .False., .False.)
437 : end do
438 : end if
439 :
440 1395 : 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 1395 : mqgridvl_old=psps%mqgrid_vl
445 : end if
446 :
447 5285 : if(ntypat_old/=dtset%ntypat.or. usepaw_old/=psps%usepaw)then
448 1392 : changed = changed + 1
449 1392 : if(idtset/=1) then
450 0 : ABI_SFREE(psps%xccc1d)
451 0 : ABI_SFREE(psps%xcctau1d)
452 : end if
453 5568 : ABI_MALLOC(psps%xccc1d,(psps%n1xccc*(1-psps%usepaw),6,dtset%ntypat))
454 4176 : ABI_MALLOC(psps%xcctau1d,(psps%n1xccc*(1-psps%usepaw),6,dtset%ntypat))
455 11824344 : psps%xccc1d = zero
456 11824344 : psps%xcctau1d = zero
457 1392 : usepaw_old=psps%usepaw
458 : end if
459 :
460 5285 : if(ntypat_old/=dtset%ntypat)then
461 1392 : changed = changed + 1
462 1392 : 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 4176 : ABI_MALLOC(psps%xcccrc,(dtset%ntypat))
468 2784 : ABI_MALLOC(psps%znucltypat,(dtset%ntypat))
469 2784 : ABI_MALLOC(psps%ziontypat,(dtset%ntypat))
470 3282 : psps%xcccrc = zero
471 3282 : psps%znucltypat = zero
472 3282 : psps%ziontypat = zero
473 1392 : ntypat_old=dtset%ntypat
474 : end if
475 :
476 12183 : psps%ziontypat(:)=dtset%ziontypat(:)
477 :
478 5285 : 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 6999 : 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 6999 : ABI_SFREE(psps%filpsp)
503 6999 : ABI_SFREE(psps%pspcod)
504 6999 : ABI_SFREE(psps%pspdat)
505 6999 : ABI_SFREE(psps%pspso)
506 6999 : ABI_SFREE(psps%pspxc)
507 6999 : ABI_SFREE(psps%title)
508 6999 : ABI_SFREE(psps%algalch)
509 6999 : ABI_SFREE(psps%mixalch)
510 6999 : ABI_SFREE(psps%ekb)
511 6999 : ABI_SFREE(psps%indlmn)
512 6999 : ABI_SFREE(psps%ffspl)
513 6999 : ABI_SFREE(psps%qgrid_ff)
514 6999 : ABI_SFREE(psps%qgrid_vl)
515 6999 : ABI_SFREE(psps%vlspl)
516 6999 : ABI_SFREE(psps%dvlspl)
517 6999 : ABI_SFREE(psps%xccc1d)
518 6999 : ABI_SFREE(psps%xcctau1d)
519 6999 : ABI_SFREE(psps%xcccrc)
520 6999 : ABI_SFREE(psps%ziontypat)
521 6999 : ABI_SFREE(psps%zionpsp)
522 6999 : ABI_SFREE(psps%znucltypat)
523 6999 : ABI_SFREE(psps%znuclpsp)
524 6999 : ABI_SFREE(psps%md5_pseudos)
525 6999 : ABI_SFREE(psps%epsatm)
526 :
527 : ! Free types.
528 6999 : call psp2params_free(psps%gth_params)
529 :
530 6999 : if (allocated(psps%nctab)) then
531 10944 : do ii=1,size(psps%nctab)
532 10944 : call nctab_free(psps%nctab(ii))
533 : end do
534 10944 : ABI_FREE(psps%nctab)
535 : end if
536 :
537 6999 : 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 3445 : 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 3445 : pspsout%dimekb = pspsin%dimekb
563 3445 : pspsout%lmnmax = pspsin%lmnmax
564 3445 : pspsout%lnmax = pspsin%lnmax
565 3445 : pspsout%mproj = pspsin%mproj
566 3445 : pspsout%mpsang = pspsin%mpsang
567 3445 : pspsout%mpspso = pspsin%mpspso
568 3445 : pspsout%mpssoang = pspsin%mpssoang
569 3445 : pspsout%mqgrid_ff = pspsin%mqgrid_ff
570 3445 : pspsout%mqgrid_vl = pspsin%mqgrid_vl
571 3445 : pspsout%mtypalch = pspsin%mtypalch
572 3445 : pspsout%npsp = pspsin%npsp
573 3445 : pspsout%npspalch = pspsin%npspalch
574 3445 : pspsout%ntypat = pspsin%ntypat
575 3445 : pspsout%ntypalch = pspsin%ntypalch
576 3445 : pspsout%ntyppure = pspsin%ntyppure
577 3445 : pspsout%n1xccc = pspsin%n1xccc
578 3445 : pspsout%optnlxccc = pspsin%optnlxccc
579 3445 : pspsout%positron = pspsin%positron
580 3445 : pspsout%usepaw = pspsin%usepaw
581 3445 : pspsout%usewvl = pspsin%usewvl
582 3445 : pspsout%useylm = pspsin%useylm
583 3445 : pspsout%nc_xccc_gspace = pspsin%nc_xccc_gspace
584 :
585 : ! logical
586 3445 : pspsout%vlspl_recipSpace = pspsin%vlspl_recipSpace
587 :
588 : ! integer allocatable
589 3445 : if (allocated(pspsin%algalch)) call alloc_copy(pspsin%algalch, pspsout%algalch)
590 3445 : if (allocated(pspsin%indlmn)) call alloc_copy(pspsin%indlmn, pspsout%indlmn)
591 3445 : if (allocated(pspsin%pspdat)) call alloc_copy(pspsin%pspdat, pspsout%pspdat)
592 3445 : if (allocated(pspsin%pspcod)) call alloc_copy(pspsin%pspcod, pspsout%pspcod)
593 3445 : if (allocated(pspsin%pspso)) call alloc_copy(pspsin%pspso, pspsout%pspso)
594 3445 : if (allocated(pspsin%pspxc)) call alloc_copy(pspsin%pspxc, pspsout%pspxc)
595 :
596 : ! real allocatable
597 3445 : if (pspsin%dimekb > 0 .and. pspsin%usepaw==0) then
598 2791 : if (allocated(pspsin%ekb)) then
599 2791 : call alloc_copy( pspsin%ekb, pspsout%ekb)
600 : end if
601 : else
602 2616 : ABI_MALLOC(pspsout%ekb,(pspsout%dimekb,pspsout%ntypat * (1 - pspsout%usepaw)))
603 654 : pspsout%ekb = zero
604 : end if
605 3445 : if (allocated(pspsin%ffspl)) call alloc_copy( pspsin%ffspl, pspsout%ffspl)
606 3445 : if (allocated(pspsin%mixalch)) call alloc_copy(pspsin%mixalch, pspsout%mixalch)
607 3445 : if (allocated(pspsin%qgrid_ff)) call alloc_copy(pspsin%qgrid_ff, pspsout%qgrid_ff)
608 3445 : if (allocated(pspsin%qgrid_vl)) call alloc_copy(pspsin%qgrid_vl, pspsout%qgrid_vl)
609 3445 : if (allocated(pspsin%vlspl)) call alloc_copy(pspsin%vlspl, pspsout%vlspl)
610 3445 : if (allocated(pspsin%dvlspl)) call alloc_copy(pspsin%dvlspl, pspsout%dvlspl)
611 :
612 3445 : if (allocated(pspsin%ziontypat)) call alloc_copy(pspsin%ziontypat, pspsout%ziontypat)
613 3445 : if (allocated(pspsin%znucltypat)) call alloc_copy(pspsin%znucltypat, pspsout%znucltypat)
614 3445 : if (allocated(pspsin%epsatm)) call alloc_copy(pspsin%epsatm,pspsout%epsatm)
615 :
616 : ! GA: Could make a check on mtypalch here
617 3445 : if (allocated(pspsin%znuclpsp)) call alloc_copy(pspsin%znuclpsp, pspsout%znuclpsp)
618 3445 : if (allocated(pspsin%zionpsp)) call alloc_copy(pspsin%zionpsp, pspsout%zionpsp)
619 :
620 3445 : if (pspsin%n1xccc > 0) then
621 2013 : if (allocated(pspsin%xcccrc)) call alloc_copy(pspsin%xcccrc, pspsout%xcccrc)
622 2013 : if (allocated(pspsin%xccc1d)) call alloc_copy(pspsin%xccc1d, pspsout%xccc1d)
623 2013 : if (allocated(pspsin%xcctau1d)) call alloc_copy(pspsin%xcctau1d, pspsout%xcctau1d)
624 : end if
625 :
626 : ! allocate and copy character strings
627 10335 : ABI_MALLOC(pspsout%filpsp,(pspsout%npsp))
628 6890 : ABI_MALLOC(pspsout%title,(pspsout%npsp))
629 10335 : ABI_MALLOC(pspsout%md5_pseudos,(pspsout%npsp))
630 7821 : do ii=1,pspsout%npsp
631 4376 : pspsout%filpsp(ii) = pspsin%filpsp(ii)
632 4376 : pspsout%title(ii) = pspsin%title(ii)
633 7821 : pspsout%md5_pseudos(ii) = pspsin%md5_pseudos(ii)
634 : end do
635 :
636 : ! allocate and copy objects
637 3445 : if (allocated(pspsin%nctab)) then
638 11740 : ABI_MALLOC(pspsout%nctab,(pspsout%ntypat))
639 2792 : if (pspsin%usepaw==0) then
640 6153 : do ii=1,pspsout%ntypat
641 6153 : call nctab_copy(pspsin%nctab(ii), pspsout%nctab(ii))
642 : end do
643 : end if
644 : end if
645 :
646 3445 : call psp2params_copy(pspsin%gth_params, pspsout%gth_params)
647 :
648 3445 : 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 2262 : 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 2262 : my_prtvol=0; if (present(prtvol)) my_prtvol=prtvol
690 2262 : mode='COLL'; if (present(mode_paral)) mode=mode_paral
691 :
692 : ! General info including spin-orbit
693 2262 : call wrtout(units,' ==== Info on pseudopotentials ==== ', mode)
694 :
695 4001 : select case (psps%usepaw)
696 : case (0)
697 1739 : 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 523 : write(msg,'(a)')
704 523 : 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 2262 : 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 2262 : ' Number of pseudopotentials .. ',psps%npsp,ch10,&
722 4524 : ' Number of types of atoms .. ',psps%ntypat
723 2262 : call wrtout(units,msg,mode)
724 :
725 2262 : if (psps%usepaw==0) then
726 3148 : select case (psps%mpspso)
727 : case (1)
728 1409 : call wrtout(units,' Scalar calculation (no spin-orbit term) ',mode)
729 : case (2)
730 : write(msg,'(3a,i3)')&
731 330 : ' Calculation with spin-orbit coupling ',ch10,&
732 660 : ' Max number of channels (spin-orbit included) ',psps%mpssoang
733 330 : call wrtout(units,msg,mode)
734 780 : do itypat=1,psps%ntypat
735 780 : if (psps%pspso(itypat) /= 1) then
736 : write(msg,'(a,i4,a,i2,a)')&
737 450 : ' - Atom type ',itypat,' has spin-orbit characteristics (pspso= ',psps%pspso(itypat),")"
738 450 : call wrtout(units,msg,mode)
739 : end if
740 : end do
741 : case default
742 1739 : ABI_ERROR(sjoin("Invalid mpspso: ", itoa(psps%mpspso)))
743 : end select
744 : else
745 1266 : select case (maxval(psps%pspso))
746 : case (0,1)
747 468 : msg=' Scalar calculation (no spin-orbit term) '
748 : case (2)
749 523 : msg=' Calculation with spin-orbit coupling '
750 : end select
751 523 : call wrtout(units,msg,mode)
752 : end if
753 :
754 : ! Info on nonlocal part
755 3943 : select case (psps%useylm)
756 : case (0)
757 1681 : msg = ' Nonlocal part applied using Legendre polynomials '
758 : case (1)
759 581 : msg = ' Nonlocal part applied using real spherical harmonics '
760 : case default
761 2262 : ABI_ERROR(sjoin("Invalid useylm: ", itoa(psps%useylm)))
762 : end select
763 2262 : call wrtout(units,msg,mode)
764 :
765 2262 : write(msg,'(a,i3)')' Max number of non-local projectors over l and type ',psps%mproj
766 2262 : call wrtout(units,msg,mode)
767 :
768 : write(msg,'(a,i3,2a,i3,2a,i3)')&
769 2262 : ' Highest angular momentum +1 ....... ',psps%mpsang,ch10,&
770 2262 : ' Max number of (l,n) components .. ',psps%lnmax, ch10,&
771 4524 : ' Max number of (l,m,n) components .. ',psps%lmnmax
772 2262 : call wrtout(units,msg,mode)
773 :
774 : !FIXME for paw n1xccc==1
775 : ! Non-linear Core correction
776 2262 : if (psps%n1xccc/=0) then
777 1367 : write(msg,'(3a,2(a,i4,a),2a)')ch10,&
778 1367 : ' Pseudo-Core Charge Info: ',ch10,&
779 1367 : ' Number of radial points for pseudo-core charge .. ',psps%n1xccc,ch10,&
780 1367 : ' XC core-correction treatment (optnlxccc) ........ ',psps%optnlxccc,ch10,&
781 2734 : ' Radius for pseudo-core charge for each type ..... ',ch10
782 1367 : call wrtout(units,msg,mode)
783 3214 : do itypat=1,psps%ntypat
784 1847 : write(msg,'(a,i4,a,f12.4)')' - Atom type ',itypat,' has pseudo-core radius .. ',psps%xcccrc(itypat)
785 3214 : call wrtout(units,msg,mode)
786 : end do
787 : end if
788 :
789 : ! Alchemical mixing
790 2262 : 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 2262 : write(msg,'(3a,a,i6,a,a,i6)')ch10,&
810 2262 : ' Info on the Q-grid used for form factors in spline form: ',ch10,&
811 2262 : ' Number of q-points for radial functions ffspl .. ',psps%mqgrid_ff,ch10,&
812 4524 : ' Number of q-points for vlspl ................... ',psps%mqgrid_vl
813 2262 : call wrtout(units,msg,mode)
814 :
815 2262 : if (psps%vlspl_recipSpace) then
816 2262 : 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 2262 : if (psps%usepaw == 0) then
821 1739 : if (psps%nc_xccc_gspace == 0) call wrtout(units,' model core charge treated in real-space', mode)
822 1739 : 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 2262 : call wrtout(units, "", mode)
827 5160 : do itypat=1,psps%ntypat
828 2898 : write(msg,'(a,i0,a,i0)')' XC functional for type ',itypat,' is ',psps%pspxc(itypat)
829 2898 : call wrtout(units,msg,mode)
830 : !write(std_out,*)psps%ziontypat(itypat),psps%znucltypat(itypat)
831 5160 : if (psps%usepaw == 0) then
832 2155 : 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 2262 : if (my_prtvol>=3) then
840 451 : do ips=1,psps%npsp
841 248 : write(std_out,*)' Pseudo number ',ips,' read from ',trim(psps%filpsp(ips))
842 248 : write(std_out,*)' Format or code ',psps%pspcod(ips)
843 248 : write(std_out,*)' Generation date ',psps%pspdat(ips)
844 451 : write(std_out,*)' Content of first line: ', trim(psps%title(ips))
845 : end do
846 : end if
847 :
848 2262 : call wrtout(units, "", mode)
849 :
850 2262 : 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 3280 : 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 3280 : real(dp), allocatable :: dummy3(:,:,:)
916 : !real(dp), allocatable :: dummy1(:)
917 : ! *************************************************************************
918 :
919 3280 : 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 3280 : with_xccc = 0
927 3280 : if (psps%n1xccc > 0) then
928 1912 : with_xccc = 1
929 1912 : with_xcctau = 1
930 : end if
931 3280 : n1xccc = max(1, psps%n1xccc)
932 :
933 3280 : if (.not. allocated(psps%xcccrc) .or. .not. allocated(psps%xccc1d) .or. psps%usepaw /= 0) then
934 1815 : with_xccc = 0
935 : end if
936 3280 : if (.not. allocated(psps%xcctau1d) .or. psps%usepaw /= 0) then
937 1815 : 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 36080 : ])
953 3280 : NCF_CHECK(ncerr)
954 :
955 : ! Define variables
956 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
957 19680 : "usepaw", "useylm", "with_xccc", "with_xcctau", "with_alch"])
958 3280 : 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 36080 : ])
973 3280 : NCF_CHECK(ncerr)
974 :
975 3280 : if (psps%usepaw == 0) then
976 2677 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("ekb", "dp", "dimekb, ntypat")))
977 : !if (with_xccc > 0) then
978 2677 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("xccc1d", "dp", "n1xccc, six, ntypat")))
979 2677 : if (with_xcctau > 0) then
980 1465 : NCF_CHECK(nctk_def_arrays(ncid, nctkarr_t("xcctau1d", "dp", "n1xccc, six, ntypat")))
981 : end if
982 2677 : 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 10708 : ])
989 2677 : NCF_CHECK(ncerr)
990 : end if
991 :
992 : ! Write data
993 3280 : NCF_CHECK(nf90_put_var(ncid, vid("ziontypat"), psps%ziontypat))
994 3280 : 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 36080 : [psps%usepaw, psps%useylm, with_xccc, with_xcctau, with_alch])
1000 3280 : NCF_CHECK(ncerr)
1001 :
1002 3280 : if (allocated(psps%pspso)) then
1003 3280 : NCF_CHECK(nf90_put_var(ncid, vid("spinorbit"), psps%pspso))
1004 : end if
1005 7443 : do ipsp=1,psps%npsp
1006 12489 : NCF_CHECK(nf90_put_var(ncid, vid("filpsp"), trim(psps%filpsp(ipsp)), start=[1, ipsp]))
1007 15769 : NCF_CHECK(nf90_put_var(ncid, vid("md5_pseudos"), trim(psps%md5_pseudos(ipsp)), start=[1, ipsp]))
1008 : end do
1009 3280 : if (allocated(psps%qgrid_vl)) then
1010 3278 : NCF_CHECK(nf90_put_var(ncid, vid("qgrid_vl"), psps%qgrid_vl))
1011 : end if
1012 3280 : if (allocated(psps%qgrid_ff)) then
1013 3278 : NCF_CHECK(nf90_put_var(ncid, vid("qgrid_ff"), psps%qgrid_ff))
1014 : end if
1015 3280 : if (allocated(psps%indlmn)) then
1016 3280 : NCF_CHECK(nf90_put_var(ncid, vid("indlmn"), psps%indlmn))
1017 : end if
1018 :
1019 : ! Local part in q-space and second derivative
1020 3280 : if (allocated(psps%vlspl)) then
1021 3278 : 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 3280 : if (allocated(psps%ffspl)) then
1027 3278 : NCF_CHECK(nf90_put_var(ncid, vid("ffspl"), psps%ffspl))
1028 : end if
1029 :
1030 3280 : if (with_xccc > 0) then
1031 :
1032 : ! Pseudo-core charge for each type of atom, on the real-space radial
1033 1465 : NCF_CHECK(nf90_put_var(ncid, vid("xcccrc"), psps%xcccrc))
1034 1465 : NCF_CHECK(nf90_put_var(ncid, vid("xccc1d"), psps%xccc1d))
1035 1465 : if (with_xcctau > 0) then
1036 1465 : 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 5957 : if (psps%usepaw == 0) then
1055 2677 : if (allocated(psps%ekb)) then
1056 2677 : NCF_CHECK(nf90_put_var(ncid, vid("ekb"), psps%ekb))
1057 : end if
1058 5890 : do itypat=1,psps%ntypat
1059 :
1060 : ! TODO Could write variables has_tvale and has_tcore
1061 3213 : if (psps%nctab(itypat)%has_tvale) then
1062 2268 : ncerr = nf90_put_var(ncid, vid("nc_tvalespl"), psps%nctab(itypat)%tvalespl, start=[1,1,itypat])
1063 567 : NCF_CHECK(ncerr)
1064 : else
1065 10584 : ABI_MALLOC(dummy3, (psps%mqgrid_vl, 2, psps%ntypat))
1066 21527348 : dummy3 = zero
1067 5292 : ncerr = nf90_put_var(ncid, vid("nc_tvalespl"), dummy3)
1068 2646 : NCF_CHECK(ncerr)
1069 2646 : ABI_FREE(dummy3)
1070 : end if
1071 5890 : if (psps%nctab(itypat)%has_tcore) then
1072 6844 : ncerr = nf90_put_var(ncid, vid("nc_tcorespl"), psps%nctab(itypat)%tcorespl, start=[1,1,itypat])
1073 1711 : NCF_CHECK(ncerr)
1074 1711 : if (with_xcctau > 0) then
1075 6844 : ncerr = nf90_put_var(ncid, vid("nc_ttaucorespl"), psps%nctab(itypat)%ttaucorespl, start=[1,1,itypat])
1076 1711 : 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 49858 : integer function vid(vname)
1094 : character(len=*),intent(in) :: vname
1095 49858 : 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 1955 : 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 5865 : ABI_MALLOC(gth_params%set,(npsp))
1322 4621 : gth_params%set(:) = .false.
1323 :
1324 : !Check array, have geometric information been filled?
1325 3910 : ABI_MALLOC(gth_params%hasGeometry,(npsp))
1326 4621 : gth_params%hasGeometry(:) = .false.
1327 :
1328 : !Coefficients for local part and projectors
1329 5865 : ABI_MALLOC(gth_params%psppar,(0:4, 0:6, npsp))
1330 116593 : gth_params%psppar = zero
1331 :
1332 : !Coefficients for spin orbit part
1333 5865 : ABI_MALLOC(gth_params%psp_k_par,(1:4, 1:3, npsp))
1334 44611 : gth_params%psp_k_par = zero
1335 :
1336 : !Different radii
1337 5865 : ABI_MALLOC(gth_params%radii_cf,(npsp, 3))
1338 15818 : gth_params%radii_cf = zero
1339 :
1340 1955 : 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 3445 : 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 3445 : if (allocated(gth_paramsin%psppar)) then
1365 3435 : call alloc_copy( gth_paramsin%psppar, gth_paramsout%psppar)
1366 : end if
1367 3445 : if (allocated(gth_paramsin%radii_cf)) then
1368 3435 : call alloc_copy( gth_paramsin%radii_cf, gth_paramsout%radii_cf)
1369 : end if
1370 3445 : if (allocated(gth_paramsin%psp_k_par)) then
1371 3435 : call alloc_copy( gth_paramsin%psp_k_par, gth_paramsout%psp_k_par)
1372 : end if
1373 3445 : if (allocated(gth_paramsin%hasGeometry)) then
1374 3435 : call alloc_copy( gth_paramsin%hasGeometry, gth_paramsout%hasGeometry)
1375 : end if
1376 3445 : if (allocated(gth_paramsin%set)) then
1377 3435 : call alloc_copy( gth_paramsin%set, gth_paramsout%set)
1378 : end if
1379 :
1380 3445 : 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 6999 : subroutine psp2params_free(gth_params)
1400 :
1401 : !Arguments ------------------------------------
1402 : class(pseudopotential_gth_type),intent(inout) :: gth_params
1403 : ! *********************************************************************
1404 :
1405 6999 : ABI_SFREE(gth_params%set)
1406 6999 : ABI_SFREE(gth_params%hasGeometry)
1407 :
1408 : ! Coefficients for local part and projectors
1409 6999 : ABI_SFREE(gth_params%psppar)
1410 :
1411 : ! Coefficients for spin orbit part
1412 6999 : ABI_SFREE(gth_params%psp_k_par)
1413 :
1414 : ! Different radii
1415 6999 : ABI_SFREE(gth_params%radii_cf)
1416 :
1417 6999 : 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 1405 : 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 1405 : 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 1405 : nctab%has_tcore = has_tcore
1448 1405 : nctab%dncdq0 = zero; nctab%d2ncdq0 = zero
1449 1405 : nctab%dtaucdq0 = zero; nctab%d2taucdq0 = zero
1450 8647733 : ABI_CALLOC(nctab%tcorespl, (mqgrid_vl, 2))
1451 8646328 : ABI_CALLOC(nctab%ttaucorespl, (mqgrid_vl, 2))
1452 :
1453 : ! tvalespl is allocated only if available.
1454 1405 : nctab%has_tvale = has_tvale
1455 1405 : nctab%dnvdq0 = zero
1456 1405 : if (has_tvale) then
1457 0 : ABI_CALLOC(nctab%tvalespl, (mqgrid_vl, 2))
1458 : end if
1459 :
1460 1405 : 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 6276 : subroutine nctab_free(nctab)
1473 :
1474 : !Arguments ------------------------------------
1475 : class(nctab_t),intent(inout) :: nctab
1476 : ! *************************************************************************
1477 :
1478 6276 : ABI_SFREE(nctab%tvalespl)
1479 6276 : ABI_SFREE(nctab%tvaletauspl)
1480 6276 : ABI_SFREE(nctab%tcorespl)
1481 6276 : ABI_SFREE(nctab%ttaucorespl)
1482 6276 : ABI_SFREE(nctab%tphi_qspl)
1483 6276 : ABI_SFREE(nctab%tphi_n)
1484 6276 : ABI_SFREE(nctab%tphi_l)
1485 6276 : ABI_SFREE(nctab%tphi_jtot)
1486 6276 : ABI_SFREE(nctab%tphi_occ)
1487 :
1488 6276 : 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 3362 : 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 3362 : nctabout%mqgrid_vl = nctabin%mqgrid_vl
1508 3362 : nctabout%has_tvale = nctabin%has_tvale
1509 3362 : nctabout%has_tcore = nctabin%has_tcore
1510 3362 : nctabout%dncdq0 = nctabin%dncdq0
1511 3362 : nctabout%d2ncdq0 = nctabin%d2ncdq0
1512 3362 : nctabout%dtaucdq0 = nctabin%dtaucdq0
1513 3362 : nctabout%d2taucdq0 = nctabin%d2taucdq0
1514 3362 : nctabout%dnvdq0 = nctabin%dnvdq0
1515 3362 : nctabout%has_tvaletau = nctabin%has_tvaletau
1516 3362 : nctabout%dnvtaudq0 = nctabin%dnvtaudq0
1517 :
1518 : ! TODO Why not check values of has_tvale and has_tcore?
1519 3362 : if (allocated(nctabin%tvalespl)) call alloc_copy(nctabin%tvalespl, nctabout%tvalespl)
1520 3362 : if (allocated(nctabin%tvaletauspl)) call alloc_copy(nctabin%tvaletauspl, nctabout%tvaletauspl)
1521 3362 : if (allocated(nctabin%tcorespl)) call alloc_copy(nctabin%tcorespl, nctabout%tcorespl)
1522 3362 : if (allocated(nctabin%ttaucorespl)) call alloc_copy(nctabin%ttaucorespl, nctabout%ttaucorespl)
1523 :
1524 3362 : 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 497 : 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 497 : nctab%has_tvale = .True.
1563 497 : if (.not. allocated(nctab%tvalespl)) then
1564 1035 : ABI_MALLOC(nctab%tvalespl, (mqgrid_vl, 2))
1565 : else
1566 152 : ABI_CHECK(size(nctab%tvalespl, dim=1) == mqgrid_vl, "wrong mqgrid_vl")
1567 : end if
1568 :
1569 497 : call pawpsp_cg(nctab%dnvdq0, d2nvdq0, mqgrid_vl, qgrid_vl, nctab%tvalespl(:,1), mesh, valr, yp1, ypn)
1570 291963 : call simp_gen(yp1, mesh%rad**2 * valr, mesh)
1571 497 : 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 497 : fact = zion / nctab%tvalespl(1,1)
1579 1510377 : nctab%tvalespl(:,1) = nctab%tvalespl(:,1) * fact
1580 :
1581 : ! Compute second derivative of tvalespl(q)
1582 497 : call paw_spline(qgrid_vl,nctab%tvalespl(:,1),mqgrid_vl,yp1,ypn,nctab%tvalespl(:,2))
1583 :
1584 497 : 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 2181 : 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 2181 : type(pawrad_type) :: core_mesh
1680 : ! *************************************************************************
1681 :
1682 2181 : ABI_CHECK(mqgrid_vl == nctab%mqgrid_vl, "wrong mqgrid_vl")
1683 :
1684 2181 : if (.not. allocated(nctab%tcorespl)) then
1685 0 : ABI_CALLOC(nctab%tcorespl, (mqgrid_vl, 2))
1686 : else
1687 2181 : ABI_CHECK(size(nctab%tcorespl, dim=1) == mqgrid_vl, "wrong mqgrid_vl")
1688 : end if
1689 :
1690 2181 : if (.not. allocated(nctab%ttaucorespl)) then
1691 0 : ABI_CALLOC(nctab%ttaucorespl, (mqgrid_vl, 2))
1692 : else
1693 2181 : 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 2181 : if (abs(xcccrc) < tol16) then
1698 1130 : nctab%has_tcore = .False.
1699 : return
1700 : end if
1701 :
1702 1051 : nctab%has_tcore = .True.
1703 : ! XCCC is given on a linear mesh.
1704 1051 : amesh = xcccrc / dble(n1xccc-1)
1705 1051 : 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 1051 : core_mesh, xccc1d(:,1), yp1, ypn)
1711 :
1712 : ! Compute second derivative of tcorespl(q)
1713 1051 : 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 1051 : core_mesh, xcctau1d(:,1), yp1, ypn)
1718 1051 : call paw_spline(qgrid_vl, nctab%ttaucorespl(:,1), mqgrid_vl, yp1, ypn, nctab%ttaucorespl(:,2))
1719 :
1720 1051 : call pawrad_free(core_mesh)
1721 :
1722 1051 : 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 : !!***
|
