-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmodule_analysis.f90
More file actions
2254 lines (1697 loc) · 78.7 KB
/
module_analysis.f90
File metadata and controls
2254 lines (1697 loc) · 78.7 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
module analysis_part
use global
use density
use fedvr3d_basis_set
use operator_spatial_space
implicit none
contains
!
! calc. the expectation of the onebody operator
!==================================================================================
!
! calc. onebody energy part < WF | onebody | WF >
!
! ______ p
! = < WF | \ { onebody }_pq * E |WF>
! / p q q
! ------
!
!= _______
! \
! \ { rho 1 }_qp * { onebody }_pq
! / p,q
! /______
!
! *
! {onebody } _pq = Integrate ( {fhi_p} onebody {fhi_q} )
!
! fhi_p is the p-th spatial orbital
!=================================================================================
subroutine onebody_expectation(onebody_spatial,rdensity1,ndim,value_return)
implicit none
integer :: ip,iq
integer,intent(in) :: ndim
complex(kind=k2),intent(in) :: onebody_spatial(ndim*(ndim+1)/2)
complex(kind=k2),intent(in) :: rdensity1(ndim,ndim)
complex(kind=k2),intent(out) :: value_return
complex(kind=k2) :: ctemp
value_return = zzero
do ip = 1, ndim
do iq =1, ndim
if(ip>=iq) then
ctemp = onebody_spatial( ia(ip) + iq )
else
ctemp = conjg(onebody_spatial( ia(iq) + ip))
endif
value_return = value_return + rdensity1(iq,ip)*ctemp
enddo
enddo
return
end subroutine onebody_expectation
!
! calc. twobody operator expectation
!
! it can be improved by symmetry
!
subroutine twobody_expectation(twobody_spatial,rdensity2,ndim,value_return)
implicit none
integer :: ip,iq,is,ir
integer,intent(in) :: ndim
complex(kind=k2),intent(in) :: twobody_spatial(ndim,ndim,ndim,ndim),&
rdensity2(ndim,ndim,ndim,ndim)
complex(kind=k2),intent(out) :: value_return
value_return = zzero
do ip =1, ndim
do iq = 1,ndim
do ir=1,ndim
do is=1,ndim
value_return = value_return + twobody_spatial(ip,ir,is,iq)*rdensity2(ip,ir,is,iq)
enddo
enddo
enddo
enddo
value_return= value_return*0.5d0
return
end subroutine twobody_expectation
!! Subroutines to perform the Fast Fourier Transform
SUBROUTINE four1(data,nn,isign)
!!Subroutine taken from Numerical recipes. It is a bit changed to adapt the code. Further information in pages 501-502 from Numerical recipes in Fortran 77 (Second edition)
INTEGER isign,nn !! isign is the sign of the exponential
REAL*8 data(1:2*nn) !! real array of length 2*nn. nn MUST be an integer power of 2
INTEGER i,istep,j,m,mmax,n
REAL*8 tempi,tempr
DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
n=2*nn
j=1
do i=1,n,2
if(j.gt.i)then
tempr=data(j)
tempi=data(j+1)
data(j)=data(i)
data(j+1)=data(i+1)
data(i)=tempr
data(i+1)=tempi
endif
m=n/2
1 if ((m.ge.2).and.(j.gt.m)) then
j=j-m
m=m/2
goto 1
endif
j=j+m
End do
mmax=2
2 if (n.gt.mmax) then
istep=2*mmax
theta=6.28318530717959d0/(isign*mmax)
wpr=-2.d0*sin(0.5d0*theta)**2
wpi=sin(theta)
wr=1.d0
wi=0.d0
do m=1,mmax,2
do i=m,n,istep
j=i+mmax
tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
data(j)=data(i)-tempr
data(j+1)=data(i+1)-tempi
data(i)=data(i)+tempr
data(i+1)=data(i+1)+tempi
continue
End do
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
End do
mmax=istep
goto 2
endif
return
END SUBROUTINE four1
!! Subroutine fastft() gives the Fast Fourier Transform for a complex function on an equidistance grid starting from t_0=0, and tstep=t_i+1-t_i.
Subroutine fastft(ft,tstep,ff,fgrid)
implicit none
!!INPUT
complex*16, intent(in) :: ft(:) !! function in time
real*8, intent(in) :: tstep !! time step
!!OUTPUT
complex*16, allocatable, intent(out) :: ff(:) !! fast fourier transform
real*8, allocatable, intent(out) :: fgrid(:) !! grid of the fast fourier transform
!!Auxiliar variables
integer :: ii, jj, kk
real*8 :: dist, fstep
real*8, allocatable :: ftaux(:) !! Auxiliar array
If (allocated(ff)) deallocate(ff)
!! To check that the number of points in grid ft is a power of two.
dist=log(dble(size(ft)))/log(2.0d0)
If (abs(dble(int(dist))-dist).gt.1d-7) then
write(*,*) 'ERROR IN fastft'
write(*,*) 'The number of temporal time step ',size(ft), ' is not a power of 2'
stop
End If
allocate(ff(1:size(ft)+1)) !! Allocate the ouput
allocate(fgrid(1:size(ft)+1)) !! Allocate the ouput
ff=cmplx(0.0d0,0.0d0) !! initialize the output
!!initialize the workarray
allocate(ftaux(1:2*size(ft)))
ftaux=0.0d0
!!rearrange the input in the workarray
Do ii=1,size(ft)
ftaux(2*ii-1)=real(ft(ii))
ftaux(2*ii)=aimag(ft(ii))
End Do
!!Fast Fourier Transform
call four1(ftaux,size(ft),-1)
!!frequency step
fstep=1.0d0/(size(ft)*tstep)
fgrid=0.0d0
!!Rearrange to the output (see pages 501-502 in Numerical recipes in Fortran)
Do ii=size(ft)+2,2*size(ft),2
jj=(ii-size(ft))/2
ff(jj)=cmplx(ftaux(ii-1),ftaux(ii))
fgrid(jj)=-1.0d0/(2.0d0*tstep)+dble(jj-1)*fstep
End Do
Do ii=2,size(ft),2
jj=(ii+size(ft))/2
ff(jj)=cmplx(ftaux(ii-1),ftaux(ii))
fgrid(jj)=-1.0d0/(2.0d0*tstep)+dble(jj-1)*fstep
End Do
ff(size(ft)+1)=cmplx(ftaux(size(ft)+1),ftaux(size(ft)+2))
fgrid(size(ft)+1)=1.0d0/(2.0d0*tstep)
!!include normalization factor
ff=ff*tstep/sqrt(2.0d0*acos(-1.0d0)) !! normalize the wavefunction
fgrid=fgrid*2.0d0*acos(-1.0d0) !! change to angular frequency
!! Delete auxiliary arrays
deallocate(ftaux)
End Subroutine fastft
!!==========================================================================
!! SUBROUTINES TO CALCULATE THE WAVEFUNCTION IN THE MOMENTUM SPACE
!!==========================================================================
!! To calculate the wavefunction in momentum space. phi is the "wavefunction times r"
!! It calculates the transformation for a single point in k radial, for all the Spherical Harmonics.
!! It is written in FE-DVR basis+Spherical Harmonics
subroutine phi_momentum(phix,rout,l,m,xglobal, weights,element,basis,fedvr_k_global,weights_k,phik)
implicit none
!! This functions calculates the function in momentum space for a radial linear momentum.
!! IMPORTANT: The input function is a radial function times a Spherical Harmonics
!! INPUT
complex(kind=k2),intent(in) :: phix(:,:)
!! phix is the radial function times the spherical harmonic
real(kind=k1), intent(in) :: rout
!! rout is the outer region
integer, intent(in), allocatable :: l(:),m(:)
!! L and M for the angular function in the argument.
real*8, intent(in), allocatable :: xglobal(:)
!! xglobal is the radial basis fedvr
real*8, intent(in), allocatable :: weights(:,:)
!! weights are the weights for the element in 1st argument and basis
!! in the element for the 2nd argument
integer, intent(in) :: element(:), basis(:)
!! element and basis in this element of the function chosen
real(kind=k1),allocatable, intent(in) :: fedvr_k_global(:)
!! all the nodes of the FE-DVR
real(kind=k1),allocatable, intent(in) :: weights_k(:,:)
!! OUTPUT
complex(kind=k2), allocatable, intent(inout) :: phik(:,:)
!! AUXILIAR VARIABLES
integer :: i, j, ij ! auxiliar variable for loops
integer :: iorbital ! auxiliar variable running for orbitals
integer :: nr !! number of radial points in position space.
integer :: nk !! number of radial points in momentum space.
integer :: nangle !! number of angular functions
integer :: norbitals !! number of orbitals
integer :: laux !! auxiliar variable for the total angular momentum
INTERFACE
!! To obtain the spherical bessel
Function besselj(n,x)
implicit none
integer :: n !! order of the Spherical Bessel Function
real*8 :: x !! argument of the Spherical Bessel Function
real*8 :: besselj !! Spherical Bessel Function
End Function besselj
END INTERFACE
nr=size(xglobal) !! number of radial points in position space.
nk=size(fedvr_k_global) !! number of radial points in momentum space.
nangle=size(phix(:,1))/size(xglobal) !! number of angular functions.
norbitals=size(phix(1,:)) !! number of orbitals
!! initialize the wavefunction in momentum space
If (allocated(phik)) then
deallocate(phik)
Else
continue
End If
allocate(phik(1:nangle*nk,1:size(phi(1,:)))) !! allocate memory for the function.
!#1st argument is the coordinate and the #2nd argument is the orbital.
phik=zzero !! initialize the momentum space function
Do iorbital=1,norbitals !! running in the orbitals
Do i=1,nangle !! Run in the spherical harmonics
laux=l(i) !! set the variable for the angular momentum.
Do j=1,nk !! Run in the radial coordinates of the momentum space.
Do ij=nr,1,-1 !! run in the radial part (a FE-DVR).
If (xglobal(ij).lt.rout) exit !! If it is in the inner region, do not take into account
!! Take into account the weights on the points.
If (ij.eq.nr) then !! If this is the last point of the grid.
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)+phix((i-1)*nr+ij,iorbital)*besselj(laux, fedvr_k_global(j)*xglobal(ij))*sqrt(weights(element(ij),basis(ij)))*xglobal(ij)
else !! this is not the last point of the grid.
If (basis(ij).lt.basis(ij+1)) then !! basis(ij) is not the last node in an element
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)+phix((i-1)*nr+ij,iorbital)*besselj(laux, fedvr_k_global(j)*xglobal(ij))*sqrt(weights(element(ij),basis(ij)))*xglobal(ij) !! We include the factor xglobal(ij) to have the reduced radial wavefuction
Else !! basis(ij) is the last node of the element.
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)+phix((i-1)*nr+ij,iorbital)*besselj(laux, fedvr_k_global(j)*xglobal(ij))*sqrt((weights(element(ij),basis(ij))+weights(element(ij+1),1)))*xglobal(ij) !! We include the factor xglobal(ij) to have the reduced radial wavefuction
Endif
End If
End Do
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt(2.0/acos(-1.0d0))*(-ci)**dble(laux)*fedvr_k_global(j) !! including the factor multiplying each element
!! Multiply by the weights to obtain the FE-DVR in the momentum space
!! Remember, the number of weights and nodes is the same than in the case of position space.
If (j.eq.nk) then !! last point of the grid
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt(weights_k(element(j),basis(j))) !! not DVR functions
else !! not the last point of the grid
If (basis(j).lt.basis(j+1)) then !! basis(j) is not the last node in an element
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt(weights_k(element(j),basis(j))) !! not DVR functions
else
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt((weights_k(element(j),basis(j))+weights_k(element(j+1),1))) !! bridge function
End If
End If
End Do
End Do
End Do
!! in phik(#1,#2) we have stored the radial reduced function in momentum space. That is, "the wavefunction times k". #1 is the FE-DVR+spherical harmonics and #2 is the orbital number.
return
End subroutine phi_momentum
!! To calculate the wavefunction in momentum space. phi is the "wavefunction times r"
!! It calculates the transformation for a single point in k radial, for all the Spherical Harmonics.
!! It is written in FE-DVR basis+Spherical Harmonics
!! Here we include a Hamming Window h(r), which multiplies the wavefunction to get rid of the tail of the wavefunction after r_out. It is defined as
!! { 1-cos(pi/2*(r-r_out)/Delta), r_out <= r <= r_out+Delta
!! h(r)= {
!! { 1, r_out+Delta< r
!! where Delta is the parameter of the Hamming window
subroutine phi_momentum_hamming(phix,rout,rend,delta,l,m,xglobal, weights,element,basis,fedvr_k_global,weights_k,phik)
implicit none
!! This functions calculates the function in momentum space for a radial linear momentum.
!! IMPORTANT: The input function is a radial function times a Spherical Harmonics
!!INPUT
complex(kind=k2),intent(in) :: phix(:,:)
!! phix is the radial function times the spherical harmonic
real(kind=k1), intent(in) :: rout
!! rout is the outer region
real(kind=k1), intent(in) :: rend
!! rend is the end of the box
real(kind=k1), intent(in) :: delta
!! delta is the parameter of the Hamming window
integer, intent(in), allocatable :: l(:),m(:)
!! L and M for the angular function in the argument.
real*8, intent(in), allocatable :: xglobal(:)
!! xglobal is the radial basis fedvr
real*8, intent(in), allocatable :: weights(:,:)
!! weights are the weights for the element in 1st argument and basis
!! in the element for the 2nd argument
integer, intent(in) :: element(:), basis(:)
!! element and basis in this element of the function chosen
real(kind=k1),allocatable, intent(in) :: fedvr_k_global(:)
!! all the nodes of the FE-DVR
real(kind=k1),allocatable, intent(in) :: weights_k(:,:)
!! OUTPUT
complex(kind=k2), allocatable, intent(inout) :: phik(:,:)
!! AUXILIAR VARIABLES
integer :: i, j, ij ! auxiliar variable for loops
integer :: iorbital ! auxiliar variable running for orbitals
integer :: nr !! number of radial points in position space.
integer :: nk !! number of radial points in momentum space.
integer :: nangle !! number of angular functions
integer :: norbitals !! number of orbitals
integer :: laux !! auxiliar variable for the total angular momentum
complex(kind=k2), allocatable :: phix_hamming(:,:)
!! phix_hamming is the radial function times the spherical harmonic times
!! the Hamming Window, h(r)
INTERFACE
!! To obtain the spherical bessel
Function besselj(n,x)
implicit none
integer :: n !! order of the Spherical Bessel Function
real*8 :: x !! argument of the Spherical Bessel Function
real*8 :: besselj !! Spherical Bessel Function
End Function besselj
END INTERFACE
nr=size(xglobal) !! number of radial points in position space.
nk=size(fedvr_k_global) !! number of radial points in momentum space.
nangle=size(phix(:,1))/size(xglobal) !! number of angular functions.
norbitals=size(phix(1,:)) !! number of orbitals
!! initialize the wavefunction in position space
allocate(phix_hamming(1:nr*nangle,1:norbitals))
phix_hamming=zzero
!! Now, we multiply the orbitals times the Hamming Window, h(r)
Do i=1,norbitals !! Running in the orbitals
Do j=1, nangle !! Running in the angular part
Do ij=1,nr !! Running in the radial part
If (xglobal(ij).gt.rout.and.xglobal(ij).lt.rout+delta) then
phix_hamming((j-1)*nr+ij,i)=(1.0d0-cos(acos(-1.0d0)*(xglobal(ij)-rout)/(2.0d0*delta)))*phix((j-1)*nr+ij,i)
Else If (xglobal(ij).gt.rend-delta.and.xglobal(ij).lt.rend) then
phix_hamming((j-1)*nr+ij,i)=cos(acos(-1.0d0)*(xglobal(ij)-rend+delta)/(2.0d0*delta))*phix((j-1)*nr+ij,i)
Else
phix_hamming((j-1)*nr+ij,i)=phix((j-1)*nr+ij,i)
End If
End Do
End Do
End Do
!! initialize the wavefunction in momentum space
If (allocated(phik)) then
deallocate(phik)
Else
continue
End If
allocate(phik(1:nangle*nk,1:size(phi(1,:)))) !! allocate memory for the function.
!#1st argument is the coordinate and the #2nd argument is the orbital.
phik=zzero !! initialize the momentum space function
Do iorbital=1,norbitals !! running in the orbitals
Do i=1,nangle !! Run in the spherical harmonics
laux=l(i) !! set the variable for the angular momentum.
Do j=1,nk !! Run in the radial coordinates of the momentum space.
Do ij=nr,1,-1 !! run in the radial part (a FE-DVR).
If (xglobal(ij).lt.rout) exit !! If it is in the inner region, do not take into account
!! Take into account the weights on the points.
If (ij.eq.nr) then !! If this is the last point of the grid.
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)+phix_hamming((i-1)*nr+ij,iorbital)*besselj(laux, fedvr_k_global(j)*xglobal(ij))*sqrt(weights(element(ij),basis(ij)))*xglobal(ij)
else !! this is not the last point of the grid.
If (basis(ij).lt.basis(ij+1)) then !! basis(ij) is not the last node in an element
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)+phix_hamming((i-1)*nr+ij,iorbital)*besselj(laux, fedvr_k_global(j)*xglobal(ij))*sqrt(weights(element(ij),basis(ij)))*xglobal(ij) !! We include the factor xglobal(ij) to have the reduced radial wavefuction
Else !! basis(ij) is the last node of the element.
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)+phix_hamming((i-1)*nr+ij,iorbital)*besselj(laux, fedvr_k_global(j)*xglobal(ij))*sqrt((weights(element(ij),basis(ij))+weights(element(ij+1),1)))*xglobal(ij) !! We include the factor xglobal(ij) to have the reduced radial wavefuction
Endif
End If
End Do
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt(2.0/acos(-1.0d0))*(-ci)**dble(laux)*fedvr_k_global(j) !! including the factor multiplying each element
!! Multiply by the weights to obtain the FE-DVR in the momentum space
!! Remember, the number of weights and nodes is the same than in the case of position space.
If (j.eq.nk) then !! last point of the grid
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt(weights_k(element(j),basis(j))) !! not DVR functions
else !! not the last point of the grid
If (basis(j).lt.basis(j+1)) then !! basis(j) is not the last node in an element
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt(weights_k(element(j),basis(j))) !! not DVR functions
else
phik((i-1)*nk+j,iorbital)=phik((i-1)*nk+j,iorbital)*sqrt((weights_k(element(j),basis(j))+weights_k(element(j+1),1))) !! bridge function
End If
End If
End Do
End Do
End Do
!! in phik(#1,#2) we have stored the radial reduced function in momentum space. That is, "the wavefunction times p". #1 is the FE-DVR+spherical harmonics and #2 is the orbital number.
deallocate(phix_hamming)
return
End subroutine phi_momentum_hamming
!! To plot the monoparticular orbitals and the radial density in momentum space in momentum space.
subroutine plot_phik_fedvr3d(phiuni,xglobal,nb_angle,element,basis,l,m, weights,rho)
implicit none
! This subroutine produces the density for the monoparticular
! functions in momentum space.
!!INPUT
complex(kind=k2), intent(in) :: phiuni(:,:)
!! phi in terms of the fedvr global basis.
real*8, intent(in), allocatable :: xglobal(:)
!! xglobal is the radial basis fedvr
real*8, intent(in), allocatable :: weights(:,:)
!! weights are the weights for the element in 1st argument and basis
!! in the element for the 2nd argument
integer, intent(in) :: element(:), basis(:)
!! element and basis in this element of the function chosen
integer, intent(in) :: nb_angle
!! number of angular functions
integer, intent(in), allocatable :: l(:),m(:)
!! L and M for the angular function in the argument.
complex(8), allocatable, intent(in) :: rho(:,:)
!! First arguments is the annihilation and the second creation
!! wavefunction
real*8, allocatable :: density(:)
!! this is the radial wavefunction squared
!!Auxiliary aspects
integer :: i,j,k,ko, counter
character*7 :: file_orbital
real*8 :: exp_k2,exp_k,exp_km1 !! <k**2>, <k>, <k**-1>
real*8 :: exp_km2
real*8 :: exp_l2, exp_m2, exp_m, norm
allocate(density(size(xglobal)))
Do i=1,size(phiuni(1,:)) !! Run in the orbitals
!Initialize the varibles
density=0.0d0 !!Density
exp_k2=0.0d0 !! <k2>
exp_k=0.0d0 !! <k>
exp_km1=0.0d0 !! <k**-1>
exp_km2=0.0d0 !! <k**-2>
exp_l2=0.0d0 !! <L**2>
exp_m2=0.0d0 !! <M**2>
exp_m=0.0d0 !! <M>
counter=0
write(file_orbital,'(a5,i2.2)') 'phik_', i
Open(122+i,file=file_orbital)
Do j=1,nb_angle !! Run for the spherical harmonics
Do k=1,size(density) !! Run for the radial points
counter=counter+1 !! Run in the global basis (radial+angular)
!! Density of the wavefunction
density(k)=density(k)+real(conjg(phiuni(counter,i))*phiuni(counter,i))
exp_l2=exp_l2+real(conjg(phiuni(counter,i))*phiuni(counter,i))*dble(l(j)*(l(j)+1))
exp_m2=exp_m2+real(conjg(phiuni(counter,i))*phiuni(counter,i))*dble(m(j)*m(j))
exp_m=exp_m+real(conjg(phiuni(counter,i))*phiuni(counter,i))*dble(m(j))
End Do
End Do
!! Calculation of the expectation values
Do k=1,size(density)
exp_k2=exp_k2+density(k)*xglobal(k)**2.0d0
exp_k=exp_k+density(k)*xglobal(k)
exp_km1=exp_km1+density(k)*xglobal(k)**(-1.0d0)
exp_km2=exp_km2+density(k)*xglobal(k)**(-2.0d0)
End Do
write(122+i,'(a29,i2.2,a3)') ' # This is the density |phik_',i,'|^2'
write(122+i,*) '# <k**2>=', exp_k2
write(122+i,*) '# <k>=', exp_k
write(122+i,*) '# <k**-1>=', exp_km1
write(122+i,*) '# <k**-2>=', exp_km2
write(122+i,*) '# <L**2>=', exp_l2
write(122+i,*) '# <M>=', exp_m
write(122+i,*) '# <M**2>=', exp_m2
write(122+i,*) ''
write(122+i,*) ''
write(122+i,*) 0.0d0,0.0d0
Do k=1,size(density)
!! Plot the wave function.
!! Take into account the weights on the points.
If (k.lt.size(density)) then
If (basis(k).lt.basis(k+1)) then !! basis(k) is not the last
!! node of the element.
write(122+i,*) xglobal(k), density(k)/weights(element(k),basis(k))
Else !! basis(k) is the last node of the element.
write(122+i,*) xglobal(k), density(k)/(weights(element(k),basis(k))+weights(element(k+1),1))
Endif
End If
End Do
Close(122+i)
End Do
!! initialize the density
density=0.0d0
!! Calculation of the one particle density.
!! It is normalized to the number of electrons.
Do i=1,size(rho(1,:)) !! loop of the orbitals
Do j=i+1,size(rho(:,1)) !! loop of the orbitals
counter=0
Do ko=1,nb_angle
Do k=1,size(density) !! run the number of points
counter=counter+1
!! Average density
density(k)=density(k)+2.0d0*real(conjg(phiuni(counter,i))*phiuni(counter,j)*rho(i,j))
End Do
End Do
End Do
!!For j=i
counter=0
Do ko=1,nb_angle
Do k=1,size(density) !! run the number of points
counter=counter+1
!! Average density
density(k)=density(k)+real(conjg(phiuni(counter,i))*phiuni(counter,i)*rho(i,i))
End Do
End Do
End Do
!! Expectation value for the density
!! Initialize the expectation values
exp_k2=0.0d0
exp_k=0.0d0
exp_km1=0.0d0
exp_km2=0.0d0
norm=0.0d0
!! Calculation of the expectation values
Do k=1,size(density)
norm=norm+density(k) !! Norm of the density
exp_k2=exp_k2+density(k)*xglobal(k)**2.0d0
exp_k=exp_k+density(k)*xglobal(k)
exp_km1=exp_km1+density(k)*xglobal(k)**(-1.0d0)
exp_km2=exp_km2+density(k)*xglobal(k)**(-2.0d0)
End Do
exp_k2=exp_k2/norm
exp_k=exp_k/norm
exp_km1=exp_km1/norm
exp_km2=exp_km2/norm
!! Opening the density file in momentum space
Open(122,file='rho_k')
!! Writting the heading of the file
write(122,*) '# <p**2>=', exp_k2
write(122,*) '# <p>=', exp_k
write(122,*) '# <p**-1>=', exp_km1
write(122,*) '# <p**-2>=', exp_km2
write(122,*) '# Norm=', norm
write(122,*) ''
write(122,*) ''
!! Prepare the density to plot it by dividing by the weights.
write(122,*) 0.0d0,0.0d0
Do k=1,size(density)
!! Plot the wave function.
!! Take into account the weights on the points.
If (k.lt.size(density)) then
If (basis(k).lt.basis(k+1)) then !! basis(k) is not the last
!! node of the element.
write(122,*) xglobal(k), density(k)/weights(element(k),basis(k))/norm
Else !! basis(k) is the last node of the element.
write(122,*) xglobal(k), density(k)/(weights(element(k),basis(k))+weights(element(k+1),1))/norm
Endif
End If
End Do
Close(122)
Deallocate(density)
End subroutine plot_phik_fedvr3d
!!====================================================================
!!
!! CALCULATION OF THE DENSITY FUNCTION INCLUDING THE ANGLES.
!!
!! The density function is written in FE-DVR + SPHERICAL HARMONICS,
!! as the orbitals, but including more spherical harmonics.
!! The coefficients include the weights of the FE-DVR, that is, to
!! perform the integrals we do not need to include them again.
!! This subroutines can be used in both, position and momentum space.
!!
!!====================================================================
!!--------------------------------------------------------------------
!! The subroutine basis_set_rho produces the appropiate basis set in
!! both FE-DVR + SPHERICAL COORDINATES
!!--------------------------------------------------------------------
Subroutine basis_set_rho3d(element,basis,lmax,mmax,label)
implicit none
!! INPUT
integer, allocatable,intent(in) :: element(:), basis(:)
!! element and basis in this element of the function chosen
integer, intent(in) :: lmax, mmax
!! maximum value of single particle angular momentum and magnetic
!! quantum number
!! OUTPUT
integer, allocatable,intent(out) :: label(:,:)
!! label(#1,1)-> FE-DVR function of #1
!! label(#1,2)-> angular momentum of #1
!! label(#1,3)-> magnetic quantum number of #1
!! AUXILIAR
integer :: i, j, ij, l, m
integer :: counter
!! First, how many elements are there in the new basis set?
counter=0
Do l=0,2*lmax !! Running in the total angular momentum
Do m=max(-l,-2*mmax),min(l,2*mmax)
Do i=1,size(element) !! FE-DVR basis set
counter=counter+1
End Do
End Do
End Do
!! Now, we can allocate the label(:,:)
allocate(label(1:counter,1:3))
label=0 !! Initialize label
!! Now, we assign the labels
counter=0
Do l=0,2*lmax !! Running in the total angular momentum
Do m=max(-l,-2*mmax),min(l,2*mmax) !! Running magnetic quantum number
Do i=1,size(element) !! FE-DVR basis set
counter=counter+1
label(counter,1)=i
label(counter,2)=l
label(counter,3)=m
End Do
End Do
End Do
!! label assigned
return
End Subroutine basis_set_rho3d
!!---------------------------------------------------------------------------
!! This subroutine calculates the 3D one density function. ***CHECK***
!!---------------------------------------------------------------------------
!!
!! Find details of the structure of the outcoming result just before
!! the subroutine basis_set_rho3d.
!!
!!---------------------------------------------------------------------------
Subroutine rho3d(phiuni,xglobal,label_phi, label_rho3d,rho,n_phi, n_all)
implicit none
!! This subroutine calculates the 3D monoparticular density and the 3D density
!! for all the electrons
!!INPUT
complex(kind=k2),allocatable, intent(in) :: phiuni(:,:)
!! phi in terms of the fedvr global basis.
real (kind=k1), allocatable,intent(in) :: xglobal(:)
!! fedvr in the radial basis
integer, allocatable, intent(in) :: label_phi(:,:) !! refered to the wavefunction
!! label_phi(#1,1)-> FE-DVR function of #1
!! label_phi(#1,2)-> angular momentum of #1
!! label_phi(#1,3)-> magnetic quantum number of #1
integer, allocatable, intent(in) :: label_rho3d(:,:) !! refered to the density
!! label_rho3d(#1,1)-> FE-DVR function of #1
!! label_rho3d(#1,2)-> angular momentum of #1
!! label_rho3d(#1,3)-> magnetic quantum number of #1
complex(8), allocatable, intent(in) :: rho(:,:)
!! First arguments is the annihilation and the second creation
!! wavefunction
!!OUTPUT
complex(kind=k2), allocatable, intent(out) :: n_phi(:,:) !! monoparticular 3D density for each orbitals
complex(kind=k2), allocatable, intent(out) :: n_all(:) !! 3D density for all the electrons
!!Auxiliary quantities
integer :: norbital !! number of orbitals.
integer :: nbasis_phi !! number of basis set functions for phi.
integer :: nbasis_rho3d !! number of basis set functions for the 3D density.
integer :: iorbital, jorbital !! Running in the orbitals
integer :: ibasis_phi, jbasis_phi, nbradial !! running in the basis of the wavefunction.
integer :: ibasis_angle, jbasis_angle, nbangle
integer :: ibasis_rho3d, jbasis_rho3d !! counters for spatial basis for the density.
integer :: lrho3d, l,lp
integer :: mrho3d, m,mp
!!Functions
INTERFACE
!! FUNCTION TO OBTAIN THE 3J SYMBOLS
!!--------------------------------------------------------
!! The function wigner3j calculates the 3J wigner symbols
!! _________________
!! / \
!! | a b c |
!! | alfa beta gamma |
!! \ _________________/
!!
!!--------------------------------------------------------
FUNCTION wigner3j(a,b,c,alfa,beta,gamma)
IMPLICIT NONE
Integer :: a, b, c
Integer :: alfa, beta, gamma
Real*8 :: wigner3j
END FUNCTION wigner3j
!! FUNCTION TO OBTAIN THE sign_m1=(-1)**m
FUNCTION sign_m1(m)
IMPLICIT NONE
Integer, INTENT(IN) :: m
REAL*8 :: sign_m1
END FUNCTION sign_m1
END INTERFACE
print*,
print*, 'Skip subroutine rho3d'
print*, 'in module_analysis.f90'
print*, 'It takes too long!!'
return
print*,
!return
norbital= size(phiuni(1,:)) !! number of orbitals
nbasis_phi=size(label_phi(:,1)) !! number of basis set functions for wavefuctions
nbasis_rho3d=size(label_rho3d(:,1)) !! number of basis set functions for density
nbangle=nbasis_phi/size(xglobal) !! number of spherical harmonics for the wavefunctions
nbradial=size(xglobal) !! number of the radial part of the FE-DVR
!! First, we allocate the array for the density of the orbitals
allocate(n_phi(1:nbasis_rho3d,1:norbital))
allocate(n_all(1:nbasis_rho3d))
!! We initialize the output
n_phi=zzero
n_all=zzero
!! Now we obtain the density functions for the orbitals.
Do ibasis_rho3d=1,nbasis_rho3d !! Running in the spatial basis for the density.
lrho3d=label_rho3d(ibasis_rho3d,2)
mrho3d=label_rho3d(ibasis_rho3d,3)
Do ibasis_angle=1,nbangle !! Run in angle for the wavefunction, since
!! the FE-DVR is the same
ibasis_phi=(ibasis_angle-1)*nbradial+label_rho3d(ibasis_rho3d,1)
!! number of basis set for the wave function of the bra.
l=label_phi(ibasis_phi,2)
m=label_phi(ibasis_phi,3)
Do jbasis_angle=1,nbangle
jbasis_phi=(jbasis_angle-1)*nbradial+label_rho3d(ibasis_rho3d,1)
!! number of basis set for the wave function of the ket.
lp=label_phi(jbasis_phi,2)
mp=label_phi(jbasis_phi,3)
!! Check the allowed couplings
If (mrho3d.ne.mp-m) cycle
!! If M\ne m-m' then cycle
If (abs(lrho3d-lp).gt.l) cycle
!! If |L-l'|>l cycle
If (lrho3d+lp.lt.l) cycle
!! If L+l'<l cycle
!! Contribution for each
Do iorbital=1, size(phiuni(1,:)) !! run or the orbitals
Do jorbital=iorbital, size(phiuni(1,:)) !! run or the orbitals
If (iorbital.eq.jorbital) then !! for the contributions coming from a_iorbital^+a_iorbital
!! Expression of the density for each state.
n_phi(ibasis_rho3d,iorbital)=n_phi(ibasis_rho3d,iorbital)+ conjg(phiuni(ibasis_phi,iorbital))*phiuni(jbasis_phi,iorbital)*sign_m1(m+mrho3d)*wigner3j(lrho3d,l,lp,0,0,0)*wigner3j(lrho3d,l,lp,-mrho3d,-m,mp)*sqrt(dble((2*lrho3d+1)*(2*l+1)*(2*lp+1))/(4.0d0*acos(-1.0d0)))
n_all(ibasis_rho3d)=n_all(ibasis_rho3d) + conjg(phiuni(ibasis_phi,iorbital))*phiuni(jbasis_phi,jorbital)*sign_m1(m+mrho3d)*wigner3j(lrho3d,l,lp,0,0,0)*wigner3j(lrho3d,l,lp,-mrho3d,-m,mp)*sqrt(dble((2*lrho3d+1)*(2*l+1)*(2*lp+1))/(4.0d0*acos(-1.0d0)))*rho(iorbital,jorbital)
Else !! for contributions coming from a_iorbital^+a_jorbital
n_all(ibasis_rho3d)=n_all(ibasis_rho3d) + 2.0d0*real(conjg(phiuni(ibasis_phi,iorbital))*phiuni(jbasis_phi,jorbital)*rho(iorbital,jorbital))*sign_m1(m+mrho3d)*wigner3j(lrho3d,l,lp,0,0,0)*wigner3j(lrho3d,l,lp,-mrho3d,-m,mp)*sqrt(dble((2*lrho3d+1)*(2*l+1)*(2*lp+1))/(4.0d0*acos(-1.0d0)))
End If
End Do
End Do
End Do
End Do
End Do
return
End Subroutine rho3d
!! Using the coupled representations subroutines
Subroutine rho3d_coupled(orb, rho,n_phi, n_all)
implicit none
!!INPUT
complex (kind=k2), allocatable, intent(in) :: orb(:,:) !! orbitals in the uncoupled basis
complex(8), allocatable, intent(in) :: rho(:,:)
!! First arguments is the annihilation and the second creation
!! wavefunction
!!OUTPUT
complex(kind=k2), allocatable, intent(out) :: n_phi(:,:) !! monoparticular 3D density for each orbitals
complex(kind=k2), allocatable, intent(out) :: n_all(:) !! 3D density for all the electrons
!!AUXILIAR VARIBALES
complex(kind=k2), allocatable :: orb_coupled(:,:,:)
!! conjg(orb(rl'm',#1))orb(rl'm',#2) -> orb_coupled(rlm,#1,#2)
integer :: i_orb, j_orb
integer :: nglobal, norbital
type (llmm_lm_type) :: llmm_lm_aux !! coefficients which link the uncouple basis with the coupled basis for Spherical Harmonics.
!! Build the coefficients to link the uncoupled basis with the
!! coupled basis with maximum accesible angular momentum.
call uncoupled_to_coupled(fedvr3d%l_max,fedvr3d%m_max,2*fedvr3d%l_max,llmm_lm_aux)
!! Calculate conjg(orb(rl'm',#1))orb(rl'm',#2)
call phikphil_to_phikl(orb, llmm_lm_aux, orb_coupled) !! we obtain the
!! orbitals in the coupled basis
!! conjg(orb(rl'm',#1))orb(rl'm',#1)->orb_coupled(rlm,#1,#2)
nglobal=size(orb_coupled(:,1,1)) !! global basis
norbital=size(orb(1,:)) !! number of orbitals
!! allocate the density associated to the orbitals