Branch data Line data Source code
1 : : !!****m* ABINIT/m_outqmc
2 : : !! NAME
3 : : !! m_outqmc
4 : : !!
5 : : !! FUNCTION
6 : : !! Interface with Cambridge quantum Monte Carlo program 'CASINO'
7 : : !! See www.tcm.phy.cam.ac.uk/~mdt26/casino.html for more details.
8 : : !! M.D.Towler (mdt26 at cam.ac.uk) November 2003
9 : : !! N.D.M.Hine (nicholas.hine at imperial.ac.uk) November 2004
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, JYR, MKV, MT, FJ)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public License, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
17 : : !!
18 : : !! SOURCE
19 : :
20 : : #if defined HAVE_CONFIG_H
21 : : #include "config.h"
22 : : #endif
23 : :
24 : : #include "abi_common.h"
25 : :
26 : : module m_outqmc
27 : :
28 : : use defs_basis
29 : : use m_errors
30 : : use m_abicore
31 : : use m_xmpi
32 : : use m_hdr
33 : : use m_dtset
34 : :
35 : :
36 : : use defs_datatypes, only : pseudopotential_type
37 : : use defs_abitypes, only : mpi_type
38 : : use m_io_tools, only : get_unit
39 : : use m_geometry, only : xred2xcart
40 : : use m_results_gs , only : results_gs_type
41 : :
42 : : implicit none
43 : :
44 : : private
45 : : public :: outqmc
46 : :
47 : : contains
48 : : !!***
49 : :
50 : : !!****f* m_outqmc/outqmc
51 : : !! NAME
52 : : !! outqmc
53 : : !!
54 : : !! FUNCTION
55 : : !! Write the wave function to a file in 'pwfn.data' format. This file can be
56 : : !! read by the Cambridge quantum Monte Carlo program 'CASINO' and used as
57 : : !! trial wave function input for a variational or diffusion Monte Carlo calculation.
58 : : !! See www.tcm.phy.cam.ac.uk/~mdt26/casino.html for more details.
59 : : !!
60 : : !! INPUTS
61 : : !! cg(2,mcg)=wavefunction coefficients
62 : : !! dtset <type(dataset_type)>=all input variables for this dataset
63 : : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
64 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space
65 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
66 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
67 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
68 : : !! mpi_enreg=information about MPI parallelization
69 : : !! npwarr(nkpt)=number of planewaves in basis and on boundary for each k
70 : : !! occ(mband*nkpt*nsppol)=occupation number for each band and k
71 : : !! psps <type(pseudopotential_type)>=all the information about psps
72 : : !! results_gs <type(results_gs_type)>=results (energy and its components,
73 : : !! forces and its components, the stress tensor) of a ground-state computation.
74 : : !!
75 : : !! OUTPUT
76 : : !! Writes the file pwfn.data
77 : : !!
78 : : !! SOURCE
79 : :
80 : 0 : subroutine outqmc(cg,dtset,eigen,gprimd,hdr,kg,mcg,mpi_enreg,npwarr,occ,psps,results_gs)
81 : :
82 : : !Arguments -------------------------------
83 : : !scalars
84 : : integer :: mcg
85 : : type(dataset_type) :: dtset
86 : : type(hdr_type) :: hdr
87 : : type(mpi_type) :: mpi_enreg
88 : : type(pseudopotential_type) :: psps
89 : : type(results_gs_type) :: results_gs
90 : : !arrays
91 : : integer :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
92 : : real(dp) :: cg(2,mcg)
93 : : real(dp) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
94 : : real(dp) :: gprimd(3,3),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
95 : :
96 : : !Local variables -------------------------
97 : : !scalars
98 : : integer,parameter :: r2s_length=80
99 : : integer :: comm,iband,iband_kpt_shift,iband_sppol_shift,icg,icg_shift,icgfull
100 : : integer :: icgfull_shift,ii,ikg,ikg_shift,ikgfull,ikpt,ikpt_shift,io
101 : : integer :: iocc,isppol,jj,me,nband1,nband2,nelecs,nkgfull,ierr
102 : : real(dp) :: norm
103 : : logical :: am_master,foundkg
104 : : character(50) :: pwfnfilename
105 : : character(500) :: message
106 : : character(80) :: dft_functional,pseudo_name,pseudo_type
107 : : character(r2s_length) :: tmpr,tmpr2,tmpr3
108 : : !arrays
109 : 0 : integer :: kgfull(3,dtset%mpw*dtset%mkmem),kgmap(dtset%mpw*dtset%mkmem)
110 : 0 : real(dp) :: gcart_qmc(3),kptscart_qmc(3,dtset%nkpt)
111 : 0 : real(dp),allocatable :: cgfull(:,:),xcart_qmc(:,:)
112 : : ! *********************************************************************
113 : :
114 : : !Go away if I am not the master node.
115 : 0 : write(message,'(a,a)')ch10,' outqmc: enter '
116 : 0 : call wrtout(ab_out,message,'PERS')
117 : :
118 : 0 : am_master=.true.
119 : : if(xmpi_paral==1 .or. mpi_enreg%paral_kgb==1)then
120 : 0 : comm=mpi_enreg%comm_cell
121 : 0 : me=xmpi_comm_rank(comm)
122 [ # # ]: 0 : if(me/=0) am_master=.false.
123 : : end if
124 : : if(.not.am_master)return
125 : :
126 [ # # ]: 0 : if(mpi_enreg%paral_spinor==1)then
127 : 0 : message = ' Parallelization over spinors is not currently supported'
128 : 0 : ABI_ERROR(message)
129 : : end if
130 : :
131 : : !write(std_out,*)ch10,'outqmc: DEBUG: dtset%ndtset = ',dtset%ndtset,ch10
132 : : !Open CASINO pwfn.data file
133 [ # # ]: 0 : if (dtset%ndtset<2)then
134 : 0 : pwfnfilename='pwfn.data'
135 : : else
136 [ # # # # ]: 0 : pwfnfilename='pwfn'//trim(i2s(dtset%jdtset))//'.data'
137 : : end if
138 [ # # ]: 0 : call wrtout(ab_out,' outqmc: will open CASINO file: '//trim(pwfnfilename),'PERS')
139 : :
140 : 0 : io = get_unit()
141 : 0 : open(io,file=pwfnfilename,form='formatted',recl=300,status='unknown',iostat=ierr)
142 : :
143 [ # # ]: 0 : if(ierr/=0)then
144 [ # # ]: 0 : ABI_ERROR("Unable to open file: "//trim(pwfnfilename))
145 : : end if
146 : :
147 : : !Check if the full set of k vectors has been used in this calculation
148 [ # # ]: 0 : if (dtset%kptopt==1) then
149 : 0 : close(io,status='delete')
150 : 0 : write(message,'(3a)')' outqmc: ERROR - kptopt=1 so k-points have been ',&
151 : 0 : & 'generated in the irreducible Brillouin Zone only. ',&
152 : 0 : & 'Set kptopt=2 to obtain full set of k-points.'
153 : 0 : ABI_ERROR(message)
154 : : end if
155 : :
156 : : !Check if the full set of G vectors has been used in this calculation
157 [ # # ]: 0 : do ikpt=1,dtset%nkpt
158 [ # # ]: 0 : if (dtset%istwfk(ikpt)/=1) then
159 : 0 : close(io,status='delete')
160 : : write(message,'(a,i5,a,i2,a,a,a,a,a)')&
161 : 0 : & ' istwfk(',ikpt,')=',dtset%istwfk(ikpt),' (ie /= 1) so some ',&
162 : 0 : & 'G-vectors may not be present.',ch10,' Set istwfk=1 for each ',&
163 : 0 : & 'k-point to obtain full set.'
164 : 0 : ABI_ERROR(message)
165 : : end if
166 : : end do !ikpt
167 : :
168 : 0 : write(message,'(a)')' outqmc: QMC trial wave function file for CASINO generated by ABINIT'
169 : 0 : call wrtout(ab_out,message,'PERS')
170 : :
171 : : !Write the required quantities to pwfn.data
172 : 0 : write(io,"('QMC trial wave function file for CASINO generated by ABINIT (www.abinit.org).')")
173 : :
174 : 0 : write(io,fmt="(/'BASIC INFO'/'----------')")
175 : 0 : write(io,fmt="('Generated by:')")
176 : 0 : write(io,fmt="(' ABINIT ',a)")trim(hdr%codvsn)
177 : 0 : write(io,fmt="('Method:'/' DFT')")
178 : :
179 : 0 : write(io,fmt="('DFT Functional:')")
180 : 0 : select case (dtset%ixc)
181 : : case(0)
182 : 0 : dft_functional='No exchange-correlation.'
183 : : case(1)
184 : 0 : dft_functional='L(S)DA (Teter/Pade parametrization)'
185 : : case(2)
186 : 0 : dft_functional='LDA (Perdew-Zunger-Ceperley-Alder parametrization)'
187 : : case(3)
188 : 0 : dft_functional='LDA (old Teter rational polynomial parametrization)'
189 : : case(4)
190 : 0 : dft_functional='LDA (Wigner)'
191 : : case(5)
192 : 0 : dft_functional='LDA (Hedin-Lundqvist)'
193 : : case(6)
194 : 0 : dft_functional='LDA (X-alpha)'
195 : : case(7)
196 : 0 : dft_functional='L(S)DA (Perdew-Wang 92)'
197 : : case(8)
198 : 0 : dft_functional='L(S)DA (Perdew-Wang 92, exchange-only)'
199 : : case(9)
200 : 0 : dft_functional='L(S)DA (Perdew-Wang 92, exchange- and RPA-correlation)'
201 : : case(10)
202 : 0 : dft_functional='Diff. between ixc=7 and 9; use with accurate RPA corr. energy'
203 : : case(11)
204 : 0 : dft_functional='GGA (Perdew-Burke-Ernzerhof)'
205 : : case(12)
206 : 0 : dft_functional='GGA (Perdew-Burke-Ernzerhof, exchange-only)'
207 : : case(13)
208 : 0 : dft_functional='GGA (potential: van Leeuwen-Baerends ; energy: Perdew-Wang 92)'
209 : : case(14)
210 : 0 : dft_functional='GGA (RPBE of Zhang and Yang)'
211 : : case(15)
212 : 0 : dft_functional='GGA (RPBE of Hammer, Hansen and Norskov)'
213 : : case(16)
214 : 0 : dft_functional='GGA (HTCH)'
215 : : case(17)
216 : 0 : dft_functional='Not defined (as of 11/2003).'
217 : : case(18)
218 : 0 : dft_functional='Not defined (as of 11/2003).'
219 : : case(19)
220 : 0 : dft_functional='Not defined (as of 11/2003).'
221 : : case(20)
222 : 0 : dft_functional='Fermi-Amaldi xc for TDDFT'
223 : : case(21)
224 : 0 : dft_functional='Fermi-Amaldi xc for TDDFT with LDA xc kernel'
225 : : case(22)
226 : 0 : dft_functional='Fermi-Amaldi xc for TDDFT with Burke-Petersilka-Gross hybrid xc kernel'
227 : : case default
228 [ # # # # : 0 : dft_functional='Unknown type.'
# # # # #
# # # # #
# # # # #
# # # #
# ]
229 : : end select
230 : 0 : write(io,"(' ABINIT ixc = ',i2,' : ',a)")dtset%ixc,trim(dft_functional)
231 : :
232 : 0 : write(io,"('Pseudopotential (of first atom type)')")
233 : 0 : select case(psps%pspcod(1))
234 : : case(1)
235 : 0 : pseudo_type='ABINIT type 1' ; pseudo_name='Troullier-Martins'
236 : : case(2)
237 : 0 : pseudo_type='ABINIT type 2' ; pseudo_name='Goedecker-Teter-Hutter (GTH)'
238 : : case(3)
239 : 0 : pseudo_type='ABINIT type 3' ; pseudo_name='Hartwigsen-Goedecker-Hutter'
240 : : case(4)
241 : 0 : pseudo_type='ABINIT type 4'
242 : 0 : pseudo_name='Teter pseudo generated using the ATOM code'
243 : : case(5)
244 : 0 : pseudo_type='ABINIT type 5'
245 : 0 : pseudo_name='"Phoney" pseudo built on a Hamman grid'
246 : : case(6)
247 : 0 : pseudo_type='ABINIT type 6'
248 : 0 : pseudo_name='Fritz-Haber Institut (Troullier Martins)'
249 : : case default
250 [ # # # # : 0 : pseudo_type='Unknown pseudopotential type (as of 11/2003).' ; pseudo_name=''
# # # ]
251 : : end select
252 [ # # ]: 0 : if(dtset%ixc<10)then
253 : 0 : write(io,"(1x,a,2x,': ',a)")trim(pseudo_type),trim(pseudo_name)
254 : : else
255 : 0 : write(io,"(1x,a,3x,': ',a)")trim(pseudo_type),trim(pseudo_name)
256 : : end if
257 : :
258 : 0 : write(io,"('Plane wave cutoff (au)')")
259 : 0 : tmpr=r2s(hdr%ecut_eff,'(f12.3)')
260 : 0 : write(io,'(1x,a)')trim(tmpr)
261 : :
262 : 0 : write(io,"('Spin polarized:')")
263 : 0 : select case(dtset%nspden)
264 : : case(1)
265 : 0 : write(io,"(' .false.')")
266 : : case(2)
267 : 0 : write(io,"(' .true.')")
268 : : case(4)
269 : 0 : close(io,status='delete')
270 : 0 : write(message,'(a)')' outqmc: ERROR - nspden=4 but CASINO cannot yet deal with non-collinear spins.'
271 : 0 : ABI_ERROR(message)
272 : : case default
273 : 0 : close(io,status='delete')
274 [ # # # # ]: 0 : ABI_ERROR('Unrecognized value of nspden.')
275 : : end select
276 : :
277 : 0 : write(io,"('Total energy (au per primitive cell)')")
278 : 0 : tmpr=r2s(results_gs%etotal,'(f24.14)')
279 : 0 : write(io,'(1x,a)')trim(tmpr)
280 : 0 : write(io,"('Kinetic energy')")
281 : 0 : tmpr=r2s(results_gs%energies%e_kinetic,'(f24.14)')
282 : 0 : write(io,'(1x,a)')trim(tmpr)
283 : 0 : write(io,"('Local potential energy (Local pseudopotential energy eei + pseudopotential core-core energy eii)')")
284 : 0 : tmpr=r2s((results_gs%energies%e_localpsp+results_gs%energies%e_corepsp),'(f24.14)')
285 : 0 : write(io,'(1x,a)')trim(tmpr)
286 : 0 : write(io,"('Non-local potential energy')")
287 : 0 : tmpr=r2s(results_gs%energies%e_nlpsp_vfock,'(f24.14)')
288 : 0 : write(io,'(1x,a)')trim(tmpr)
289 : 0 : write(io,"('Electron-electron energy (Hartree Energy + Exchange-Correlation Energy)')")
290 : 0 : tmpr=r2s((results_gs%energies%e_hartree+results_gs%energies%e_xc),'(f24.14)')
291 : 0 : write(io,'(1x,a)')trim(tmpr)
292 : 0 : write(io,"('Ion-ion energy')")
293 : 0 : tmpr=r2s(results_gs%energies%e_ewald,'(f24.14)')
294 : 0 : write(io,'(1x,a)')trim(tmpr)
295 : 0 : write(io,"('Number of electrons per primitive cell')")
296 : :
297 : 0 : nelecs=0
298 [ # # ]: 0 : do ii=1,dtset%natom
299 : 0 : nelecs=nelecs+psps%ziontypat(dtset%typat(ii))
300 : : end do
301 : :
302 : 0 : write(io,'(1x,i3)')nelecs
303 : 0 : write(io,*)
304 : 0 : write(io,"('GEOMETRY'/'--------')")
305 : 0 : write(io,"('Number of atoms per primitive cell')")
306 : 0 : write(io,'(1x,i3)')dtset%natom
307 : 0 : write(io,"('Atomic numbers and positions of atoms (au)')")
308 : :
309 [ # # # # : 0 : ABI_MALLOC(xcart_qmc,(3,dtset%natom))
# # # # #
# # # ]
310 : 0 : call xred2xcart(dtset%natom,hdr%rprimd,xcart_qmc,hdr%xred)
311 [ # # ]: 0 : do ii=1,dtset%natom
312 : 0 : tmpr=r2s(xcart_qmc(1,ii),'(f24.14)')
313 : 0 : tmpr2=r2s(xcart_qmc(2,ii),'(f24.14)')
314 : 0 : tmpr3=r2s(xcart_qmc(3,ii),'(f24.14)')
315 : 0 : jj=psps%znucltypat(dtset%typat(ii))
316 : 0 : write(io,'(1x,i2,3(1x,a))')jj,trim(tmpr),trim(tmpr2),trim(tmpr3)
317 : : end do
318 [ # # ]: 0 : ABI_FREE(xcart_qmc)
319 : :
320 : 0 : write(io,"('Primitive lattice vectors (au)')")
321 [ # # ]: 0 : do ii=1,3
322 : 0 : tmpr=r2s(hdr%rprimd(1,ii),'(f24.14)')
323 : 0 : tmpr2=r2s(hdr%rprimd(2,ii),'(f24.14)')
324 : 0 : tmpr3=r2s(hdr%rprimd(3,ii),'(f24.14)')
325 : 0 : write(io,'(3(1x,a))')trim(tmpr),trim(tmpr2),trim(tmpr3)
326 : : end do
327 : :
328 : : !Copy the G vectors for the first k point into kgfull
329 : 0 : ikgfull=0
330 [ # # ]: 0 : do ikg=1,npwarr(1)
331 : 0 : ikgfull=ikgfull+1
332 [ # # ]: 0 : kgfull(1:3,ikgfull) = kg(1:3,ikg)
333 : 0 : kgmap(ikg)=ikgfull
334 : : end do
335 : 0 : ikg_shift = npwarr(1)
336 : : !Go through the other k points and look for any G vectors that haven't
337 : : !already been found and add them to the end of kgfull
338 [ # # ]: 0 : do ikpt=2,dtset%nkpt
339 [ # # ]: 0 : do ikg=ikg_shift,ikg_shift+npwarr(ikpt)
340 : 0 : foundkg = .false.
341 [ # # ]: 0 : do ii=1,ikgfull
342 : : if(kg(1,ikg)==kgfull(1,ii).and.kg(2,ikg)==kgfull(2,ii) &
343 [ # # # # : 0 : & .and.kg(3,ikg)==kgfull(3,ii)) then
# # ]
344 : 0 : foundkg=.true.
345 : 0 : kgmap(ikg)=ii
346 : : exit
347 : : end if
348 : : end do
349 : 0 : if(.not.foundkg)then
350 : 0 : ikgfull=ikgfull+1
351 [ # # ]: 0 : kgfull(1:3,ikgfull)=kg(1:3,ikg)
352 : 0 : kgmap(ikg)=ikgfull
353 : : end if
354 : : end do
355 : 0 : ikg_shift=ikg_shift+npwarr(ikpt)
356 : : end do
357 : 0 : nkgfull=ikgfull
358 : :
359 : 0 : write(io,*)
360 : 0 : write(io,"('G VECTORS'/'---------')")
361 : 0 : write(io,"('Number of G-vectors')")
362 : 0 : write(io,'(1x,i7)')nkgfull
363 : 0 : write(io,"('Gx Gy Gz (au)')")
364 : :
365 [ # # ]: 0 : do ikgfull=1,nkgfull
366 : : gcart_qmc=2*pi*(kgfull(1,ikgfull)*gprimd(1:3,1)&
367 [ # # ]: 0 : & +kgfull(2,ikgfull)*gprimd(1:3,2)+kgfull(3,ikgfull)*gprimd(1:3,3))
368 : 0 : write(io,*)gcart_qmc(1:3) ! '(3e26.16)'
369 : : end do
370 : :
371 : : !Populate the cgfull array, using the mapping in kgmap between the
372 : : !coefficients for kg in the per-kpoint list and the ones in the full list
373 : : !The number of xxx_shift's might seem excessive but the re-ordering of the
374 : : !list from (spin, kpt, band, kg) to (kpt, spin, band, kgfull) is quite
375 : : !complicated
376 [ # # # # : 0 : ABI_MALLOC(cgfull,(2,nkgfull*dtset%nspinor*dtset%nsppol*dtset%mband*dtset%nkpt))
# # # # #
# # # ]
377 [ # # # # ]: 0 : cgfull(1:2,1:nkgfull*dtset%nspinor*dtset%nsppol*dtset%mband*dtset%nkpt)=0
378 : 0 : icg_shift=1
379 [ # # ]: 0 : do isppol=1,dtset%nsppol
380 : 0 : ikg_shift=1
381 [ # # ]: 0 : if(isppol==2)then
382 : : ! Go back to the beginning of cgfull but skip the first set of isppol=1 bands
383 : 0 : icgfull_shift=nkgfull*dtset%nband(1)
384 : 0 : ikpt_shift=dtset%nkpt
385 : : else
386 : : icgfull_shift=0 ! Start at the beginning of cgfull
387 : : ikpt_shift=0
388 : : end if
389 [ # # ]: 0 : do ikpt=1,dtset%nkpt
390 [ # # ]: 0 : do iband=1,dtset%nband(ikpt+ikpt_shift)
391 : 0 : ikg=ikg_shift
392 : : ! Find the index in Abinit's coefficient list
393 [ # # ]: 0 : do icg=icg_shift,icg_shift+npwarr(ikpt)-1
394 : : ! Map it to an index in the full CASINO list with the mapping recorded
395 : : ! when kgfull was read in
396 : 0 : icgfull = kgmap(ikg)+icgfull_shift
397 [ # # ]: 0 : cgfull(1:2,icgfull)=cg(1:2,icg)
398 : 0 : ikg=ikg+1
399 : : end do !icg
400 : 0 : icg_shift=icg_shift+npwarr(ikpt)
401 : 0 : icgfull_shift=icgfull_shift+nkgfull
402 : : end do !iband
403 [ # # ]: 0 : if(isppol==2)then
404 : : ! Skip the isppol==1 bands
405 : 0 : icgfull_shift=icgfull_shift+nkgfull*dtset%nband(ikpt)
406 : : else
407 [ # # ]: 0 : if(dtset%nsppol==2)then
408 : : ! Skip the isppol==2 bands
409 : 0 : icgfull_shift=icgfull_shift+nkgfull*dtset%nband(ikpt+dtset%nkpt)
410 : : else
411 : : icgfull_shift=icgfull_shift ! Nothing to be skipped
412 : : end if
413 : : end if
414 : 0 : ikg_shift=ikg_shift+npwarr(ikpt)
415 : : end do !ikpt
416 : : end do !isppol
417 : :
418 : : !See if each orbital is normalised and check for integer occupancy of orbitals.
419 : : !These are checked by CASINO and it will complain if they are not as expected.
420 : 0 : icgfull_shift=1
421 : 0 : ii=0
422 : 0 : iocc=1
423 : :
424 [ # # ]: 0 : do ikpt=1,dtset%nkpt
425 [ # # ]: 0 : do isppol=1,dtset%nsppol
426 [ # # ]: 0 : do iband=1,dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
427 [ # # ]: 0 : if(occ(iocc)/=int(occ(iocc)))then
428 : : write(message,'(a,i5,a,i1,a,i5,a,f11.8,a)')&
429 : 0 : & 'Non-integer occupation number for kpt ',ikpt,', sppol ',isppol,', band ',iband,': occ=',occ(iocc),'.'
430 : 0 : ABI_WARNING(message)
431 : : end if
432 : 0 : iocc=iocc+1
433 : 0 : norm=0
434 [ # # ]: 0 : do icgfull=icgfull_shift,icgfull_shift+nkgfull-1
435 : 0 : norm=norm+cgfull(1,icgfull)**2+cgfull(2,icgfull)**2
436 : : end do !icgfull
437 : 0 : icgfull_shift=icgfull_shift+nkgfull
438 [ # # # # ]: 0 : if((norm<0.999).or.(norm>1.001))then
439 : : write(message,'(a,i5,a,i1,a,i5,a,f11.8,a)')&
440 : 0 : & 'Incorrectly normalised orbital for kpt ',ikpt,', sppol ',isppol,', band ',iband,': norm=',norm,'.'
441 : 0 : ABI_WARNING(message)
442 : : end if
443 : : end do !iband
444 : : end do !isppol
445 : : end do !ikpt
446 : :
447 : 0 : write(io,*)
448 : 0 : write(io,"('WAVE FUNCTION'/'-------------')")
449 : 0 : write(io,"('Number of k-points')")
450 : 0 : write(io,'(1x,i5)') dtset%nkpt
451 : :
452 [ # # ]: 0 : do ikpt=1,dtset%nkpt
453 : : kptscart_qmc(1:3,ikpt)=2*pi*(dtset%kptns(1,ikpt)*gprimd(1:3,1)&
454 [ # # ]: 0 : & +dtset%kptns(2,ikpt)*gprimd(1:3,2)+dtset%kptns(3,ikpt)*gprimd(1:3,3))
455 : : end do !ikpt
456 : :
457 : 0 : iband_kpt_shift=0
458 : 0 : icg_shift=0
459 : 0 : icgfull_shift=0
460 [ # # ]: 0 : do ikpt=1,dtset%nkpt
461 : 0 : write(io,"('k-point # ; # of bands (up spin/down spin) ; k-point coords (au)')")
462 [ # # ]: 0 : if(dtset%nsppol==2)then
463 : 0 : nband1=dtset%nband(ikpt)
464 : 0 : nband2=dtset%nband(ikpt+dtset%nkpt)
465 : : else
466 : 0 : nband1=dtset%nband(ikpt)
467 : 0 : nband2=0
468 : : end if
469 : 0 : tmpr=r2s(kptscart_qmc(1,ikpt),'(f24.14)')
470 : 0 : tmpr2=r2s(kptscart_qmc(2,ikpt),'(f24.14)')
471 : 0 : tmpr3=r2s(kptscart_qmc(3,ikpt),'(f24.14)')
472 : 0 : write(io,'(3(1x,i5),3(1x,a))')ikpt,nband1,nband2,trim(tmpr),trim(tmpr2), &
473 : 0 : & trim(tmpr3)
474 [ # # ]: 0 : do isppol=1,dtset%nsppol
475 [ # # ]: 0 : if (isppol==2) then
476 [ # # ]: 0 : iband_sppol_shift=sum(dtset%nband(1:dtset%nkpt))
477 [ # # ]: 0 : iband_kpt_shift=sum(dtset%nband((dtset%nkpt+1):(dtset%nkpt+ikpt-1)))
478 : : else
479 : 0 : iband_sppol_shift=0
480 [ # # ]: 0 : iband_kpt_shift=sum(dtset%nband(1:(ikpt-1)))
481 : : end if
482 [ # # ]: 0 : do iband=1,dtset%nband(ikpt)
483 : 0 : write(io,"('Band, spin, eigenvalue (au)')")
484 : 0 : tmpr=r2s(eigen(iband_kpt_shift+iband_sppol_shift+iband),'(f24.14)')
485 : 0 : write(io,'(2(1x,i5),1x,a)')iband,isppol,tmpr
486 : 0 : write(io,"('Eigenvector coefficients')")
487 [ # # ]: 0 : do icgfull=1,nkgfull
488 : 0 : write(io,"(1x,'(',e23.16,',',e23.16,')')")cgfull(1:2,icgfull+icgfull_shift)
489 : : end do !icgfull
490 : 0 : icgfull_shift=icgfull_shift+nkgfull
491 : : end do !iband
492 : : end do !isppol
493 : : end do !ikpt
494 : :
495 : 0 : close(io)
496 : :
497 : 0 : write(message,'(a,a)')' outqmc: done with writing of QMC trial wave function file for CASINO',ch10
498 : 0 : call wrtout(ab_out,message,'PERS')
499 : :
500 [ # # ]: 0 : end subroutine outqmc
501 : : !!***
502 : :
503 : : !!****f* m_outqmc/i2s
504 : : !! NAME
505 : : !! i2s
506 : : !!
507 : : !! FUNCTION
508 : : !! Convert integers to left justified strings that can be printed in the
509 : : !! middle of a sentence without introducing large amounts of white space.
510 : : !! Calling routine is intended to include something like:
511 : : !! integer i
512 : : !! i=12
513 : : !! write(std_out,*)'Integer number ',trim(i2s(i)),' with words at the end.'
514 : : !!
515 : : !! INPUTS
516 : : !!
517 : : !! OUTPUT
518 : : !!
519 : : !! SOURCE
520 : :
521 : 0 : function i2s(n)
522 : :
523 : : !Arguments ----------------------
524 : : integer, intent(in) :: n
525 : : character(len=20) :: i2s
526 : :
527 : : !Local variables ----------------
528 : : integer :: i,j
529 : : character :: tmp,sign
530 : :
531 : : ! *********************************************************************
532 : :
533 [ # # ]: 0 : if(n==0)then
534 : 0 : i2s='0' ; return
535 : : end if
536 [ # # ]: 0 : sign=' ' ; if(n<0)sign='-'
537 : :
538 [ # # ]: 0 : do i=1,len(i2s)
539 : 0 : i2s(i:i)=' '
540 : : end do
541 : :
542 : 0 : i=abs(n)
543 [ # # ]: 0 : do j=1,len(i2s)
544 [ # # ]: 0 : if(i==0)exit
545 : 0 : i2s(j:j)=achar(ichar('0')+mod(i,10))
546 : 0 : i=i/10
547 : : end do
548 : :
549 : 0 : i=1 ; j=len_trim(i2s)
550 : 0 : do
551 [ # # ]: 0 : if(i>=j)exit
552 : 0 : tmp=i2s(j:j)
553 : 0 : i2s(j:j)=i2s(i:i)
554 : 0 : i2s(i:i)=tmp
555 : 0 : i=i+1
556 : 0 : j=j-1
557 : : end do
558 : :
559 : 0 : i2s=trim(sign)//i2s
560 : :
561 : : end function i2s
562 : : !!***
563 : :
564 : : !!****f* ABINIT/r2s
565 : : !! NAME
566 : : !! r2s
567 : : !!
568 : : !! FUNCTION
569 : : !! Converts real variable with arbitrary format to string that can be
570 : : !! trimmed and printed in the middle of a sentence without introducing
571 : : !! large amounts of white space, as you would if you did
572 : : !! write(std_out,'(f12.6)')12.0 or similar. Note you need to pass through the
573 : : !! format string e.g. f12.6.
574 : : !!
575 : : !! Calling routine is intended to include something like:
576 : : !! USE utilities
577 : : !! REAL(dp) r
578 : : !! r=12._dp
579 : : !! tmpr=r2s(r,'(f12.6)')
580 : : !! write(std_out,*)'Real number ',trim(tmpr),' with words at the end.'
581 : : !!
582 : : !! Note : DON'T USE R2S IN A WRITE STATEMENT SINCE THIS IS ILLEGAL
583 : : !! IN FORTRAN90 (ALTHOUGH NOT IN FORTRAN200X). IF ANYONE HAS TIME, FEEL
584 : : !! FREE TO WRITE A VERSION OF THIS WHICH ISN'T ILLEGAL - SIMILAR TO
585 : : !! I2S ABOVE - SO THAT PEOPLE WHO HAVEN'T READ THIS NOTE DON'T FEEL
586 : : !! TEMPTED TO CALL R2S IN A WRITE STATEMENT.
587 : : !!
588 : : !! INPUTS
589 : : !!
590 : : !! OUTPUT
591 : : !!
592 : : !! SIDE EFFECTS
593 : : !!
594 : : !! SOURCE
595 : :
596 : 0 : function r2s(r,real_format)
597 : :
598 : : !Arguments -------------------------
599 : : real(dp),intent(in) :: r
600 : : character(len=*),intent(in) :: real_format
601 : : character(len=80) :: r2s
602 : :
603 : : ! *********************************************************************
604 : :
605 : : if(len(r2s)>0)then
606 : 0 : write(r2s,real_format)r
607 : 0 : r2s=adjustl(r2s)
608 : : end if
609 : :
610 : 0 : end function r2s
611 : : !!***
612 : :
613 : : end module m_outqmc
614 : : !!***
|
