-
Notifications
You must be signed in to change notification settings - Fork 5
/
extpot.F
1186 lines (1061 loc) · 39.7 KB
/
extpot.F
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
#include "symbol.inc"
!****************** module use external potential *********************
! Humbly written by Kuang... Apr-2014
! This whole code is only tested with PAW
!***********************************************************************
module mextpot
use prec
implicit none
type extpot
logical :: lextpot, dextpot
integer :: nx, ny, nz ! should be equal to the dimension of potential grid
real(q) :: FCORR
endtype extpot
type (extpot) EXTPT
RGRID, allocatable :: EXTPT_GRID(:,:,:) ! the external potential grid in normal layout
RGRID, allocatable :: VLM(:,:,:) ! external potential for each channel of each ion
RGRID, allocatable :: SCX(:), SCY(:), SCZ(:)
RGRID, allocatable :: SGRID(:,:,:) ! the B-spline interpolation coefficients
! parameters for real space VLM projection
integer, parameter :: KMAX=20, di_smear=5, NGRID=50
real(q), parameter :: Delta = 0.1_q
contains
function EXTPT_LEXTPOT()
logical :: EXTPT_LEXTPOT
EXTPT_LEXTPOT=EXTPT%lextpot
end function EXTPT_LEXTPOT
function EXTPT_DEXTPOT()
logical :: EXTPT_DEXTPOT
EXTPT_DEXTPOT=EXTPT%dextpot
end function EXTPT_DEXTPOT
! read INCAR FILE for external potential information
subroutine EXTPT_READER(NGXF,NGYF,NGZF,IU0,IU5,IU6)
use base
use reader_tags
use vaspxml
implicit none
integer :: NGXF, NGYF, NGZF, IU0, IU5, IU6, IU_POT=99
! local variables
integer :: ix, iy, iz, ntot
logical :: lopen, LDUM
INTEGER IDUM, IERR, N
REAL(q) RDUM
COMPLEX(q) CDUM
CHARACTER(1) :: CHARAC
lopen = .false.
OPEN(UNIT=IU5,FILE=INCAR,STATUS='OLD')
EXTPT%lextpot = .false.
EXTPT%dextpot = .false. ! initialiting logical dextpot
EXTPT%FCORR = 0.67_q
EXTPT%nx = -1
EXTPT%ny = -1
EXTPT%nz = -1
CALL PROCESS_INCAR(INCAR_F, 'LEXTPOT', EXTPT%lextpot)
CALL PROCESS_INCAR(INCAR_F, 'DEXTPOT', EXTPT%dextpot)
if( EXTPT%lextpot ) then
CALL PROCESS_INCAR(INCAR_F, 'FCORR', EXTPT%FCORR)
endif
if( EXTPT%lextpot ) then
open(unit=IU_POT,file='EXTPOT',status='old')
read(IU_POT,*) EXTPT%nx, EXTPT%ny, EXTPT%nz
! check grid size consistency
if ( EXTPT%nx /= NGXF .or. EXTPT%ny /= NGYF .or. EXTPT%nz /= NGZF ) then
write(IU6,*) 'ERROR: external potential grid size mismatch:'
write(IU6,*) ' NGXF, NGYF, NGZF:', NGXF, NGYF, NGZF
write(IU6,*) ' NX, NY, NZ:', EXTPT%nx, EXTPT%ny, EXTPT%nx
write(IU0,*) 'ERROR: external potential grid size mismatch:'
write(IU0,*) ' NGXF, NGYF, NGZF:', NGXF, NGYF, NGZF
write(IU0,*) ' NX, NY, NZ:', EXTPT%nx, EXTPT%ny, EXTPT%nx
stop
endif
allocate(EXTPT_GRID(EXTPT%nx, EXTPT%ny, EXTPT%nz))
! follow the tradition of CHGCAR format, ix be the fastest index
read(IU_POT,*) (((EXTPT_GRID(ix,iy,iz), ix=1,EXTPT%nx), &
iy=1,EXTPT%ny), iz=1,EXTPT%nz)
close(IU_POT)
endif
return
end subroutine EXTPT_READER
!******************** SUBROUTINE EXTERNAL_POT_ADD **********************
! A wrapper deal with the wierd variable type compatibility issue...
! In realmode, CVTOT allocated as complex array, but used as real grid
!***********************************************************************
SUBROUTINE EXTPT_EXTERNAL_POT_ADD(GRIDC, LATT_CUR, CVTOT)
USE prec
USE base
USE lattice
USE mpimy
USE mgrid
USE poscar
USE constant
IMPLICIT COMPLEX(q) (C)
IMPLICIT REAL(q) (A-B,D-H,O-Z)
TYPE (grid_3d) GRIDC
TYPE (latt) LATT_CUR
COMPLEX(q) :: CVTOT(GRIDC%MPLWV)
call EXTERNAL_POT_ADD_(GRIDC, LATT_CUR, CVTOT(1))
end subroutine EXTPT_EXTERNAL_POT_ADD
subroutine EXTPT_EXTERNAL_POT_ADD_PAW(POTAE, POT, NDIM, NCDIJ, LMMAX, NIP)
USE prec
USE base
USE lattice
USE mpimy
USE mgrid
USE poscar
USE constant
implicit none
INTEGER NDIM, NCDIJ, LMMAX, NIP
REAL(q) POT(:,:,:), POTAE(:,:,:)
integer ISP, LM, II
do ISP = 1, NCDIJ
do LM = 1, LMMAX
do II = 1, NDIM
POTAE(II,LM,ISP) = POTAE(II,LM,ISP) + VLM(II,LM,NIP) !+1.0_q
POT(II,LM,ISP) = POT(II,LM,ISP) + VLM(II,LM,NIP) !+1.0_q
enddo
enddo
enddo
end subroutine EXTPT_EXTERNAL_POT_ADD_PAW
!*********************** SUBROUTINE CALC_VLM ***************************
! This subroutine calculates the LM component of the external potential
! This code was written assuming realmode and NGZhalf = .True.
!***********************************************************************
SUBROUTINE EXTPT_CALC_VLM(LATT_CUR, GRIDC, T_INFO, P, LMDIM, WDES, IU0)
use prec
use base
use lattice
use mpimy
use constant
use mgrid
use poscar
use pseudo
use paw
use wave
use us
use asa
implicit none
type (grid_3d) GRIDC
type (latt) LATT_CUR
type (type_info) T_INFO
type (potcar),TARGET:: P(T_INFO%NTYP)
type (wavedes) WDES
integer :: LMDIM, IU0
! local variables
type (potcar),POINTER:: PP
complex(q) :: PHASE
real(q), allocatable :: CVEXT(:)
real(q) :: RCUT, R, RI(3),RVEC(3)
real(q) :: CTMP, RMAX, RSWITCH, val, RPREV, CTMPPREV
integer, parameter :: nprint=8
real(q) :: VLMK(KMAX), dV
integer :: NXMAX, NYMAX, NZMAX, NG, ofile=666
integer :: NI, NDIM, NT, NIP, LMMAX, LYMAX, L, M, LM, K, RNMAX, NR
integer :: iprint, IRDMAX, INDMAX, IND
integer, external :: MAXL_AUG, MAXL1
RGRID, allocatable :: VLMFULL(:,:,:)
integer, allocatable :: NLI(:)
real(q),allocatable :: YLM(:,:), DIST(:), XS(:), YS(:), ZS(:)
allocate(CVEXT(DIMREAL(GRIDC%MPLWV)))
NXMAX = EXTPT%nx
NYMAX = EXTPT%ny
NZMAX = EXTPT%nz
! copy the external potential to CVEXT, change normal layout to VASP layout
CVEXT = 0.0_q
call EXTERNAL_POT_ADD_(GRIDC, LATT_CUR, CVEXT(1))
! set up the maximum radial grid dimension, copied from SET_DD_PAW
NDIM=0
DO NT=1,T_INFO%NTYP
IF (ASSOCIATED(P(NT)%QPAW)) THEN
NDIM=MAX(NDIM, P(NT)%R%NMAX)
END IF
ENDDO
IF (NDIM == 0) RETURN
LYMAX =MAXL_AUG(T_INFO%NTYP,P)
LMMAX=(LYMAX+1)**2
! allocate the VLM array, distributed over COMM_INB
if (.not. allocated(VLM)) allocate(VLM(NDIM, LMMAX, WDES%NIONS))
VLM = 0.0_q
! allocate the full matrix
allocate(VLMFULL(NDIM,LMMAX,T_INFO%NIONS))
! Initialize B-spline interpolation
allocate( SGRID(NXMAX,NYMAX,NZMAX) )
allocate(SCX(-di_smear:di_smear), SCY(-di_smear:di_smear), SCZ(-di_smear:di_smear))
! 1D interpolate for delta grid, used later as basis for 3D B-spline
call CALC_Sdelta_COEFF(NXMAX, SCX, di_smear)
call CALC_Sdelta_COEFF(NYMAX, SCY, di_smear)
call CALC_Sdelta_COEFF(NZMAX, SCZ, di_smear)
call BSPLINE_INTERPOLATE_COEFF( EXTPT_GRID, SGRID, NXMAX, NYMAX, NZMAX, di_smear, &
SCX, SCY, SCZ, GRIDC )
! allocate temporary arrays
#ifdef MPI
IRDMAX = NGRID*NGRID*NGRID/GRIDC%COMM%NCPU+1
#else
IRDMAX = NGRID*NGRID*NGRID+1
#endif
allocate( YLM(IRDMAX,LMMAX), DIST(IRDMAX), XS(IRDMAX), &
& YS(IRDMAX), ZS(IRDMAX) )
! set up YLM for the fine integration grid, this is shared by all atom
! centers
call SETYLM_FINE_GRID( GRIDC, NGRID, YLM, LYMAX, LMMAX, XS, YS, ZS, DIST &
, IRDMAX, INDMAX, .True.)
! start projection
VLMFULL = 0.0_q
! do ion loops
ion1: do NI=1,T_INFO%NIONS
RI(:) = T_INFO%POSION(1:3,NI) ! ion position in direct coordinate
NT=T_INFO%ITYP(NI)
PP=> PP_POINTER(P, NI, NT)
RNMAX =PP%R%NMAX
RMAX = PP%R%R(RNMAX)
RCUT = RMAX + Delta ! Want to use RMAX+D as real cutoff
dV = (RCUT*2/NGRID)**3
LYMAX =MAXL1(PP)*2
! for each VLM
do L = 0, LYMAX
do M = 0, L*2
LM=L*L+M+1
VLMK(:) = 0.0_q
do K = 1, KMAX ! radial expansion
do IND = 1, INDMAX
! \sum YLM *RK *V(r) *dV
! Note this summation is distributed over GRIDC%COMM
R = DIST(IND)*RCUT ! in Angstrom
if ( R <= RCUT ) then
RVEC(1) = XS(IND)*R
RVEC(2) = YS(IND)*R
RVEC(3) = ZS(IND)*R
RVEC = matmul( transpose(LATT_CUR%B(1:3,1:3)), RVEC(1:3) )
RVEC(:) = RVEC(:) + RI(:) ! grid coordinate in direct
val = bspline( RVEC, SGRID, NXMAX, NYMAX, NZMAX )
VLMK(K) = VLMK(K) + dV* FSWITCH(R,RMAX,Delta)* &
& YLM(IND,LM)*RK(R,K,RCUT)*val
endif
enddo
enddo
do NR = RNMAX, 1, -1
R = PP%R%R(NR)
if ( R>RMAX/KMAX .or. L==0 ) then
CTMP = 0.0_q
do K = 1, KMAX
CTMP = CTMP + RK(R,K,RCUT) * VLMK(K)
enddo
CTMPPREV = CTMP
RPREV = R
else
CTMP = CTMPPREV * (R/RPREV)**L
endif
VLMFULL(NR,LM,NI) = CTMP
enddo
enddo
enddo
enddo ion1
! Reduce VLM over GRIDC%COMM
CALLMPI( M_sum_d(GRIDC%COMM, VLMFULL, NDIM*LMMAX*T_INFO%NIONS) )
! distribute the full matrix over COMM_INB to save memory space
VLM = 0.0_q
ion2: do NI=1,T_INFO%NIONS
! local storage index
! if this element is not treated locally CYCLE
IF (.NOT. DO_LOCAL(NI)) cycle ion2
NIP=NI_LOCAL(NI, WDES%COMM_INB)
do L = 0, LYMAX
do M = 0, L*2
LM=L*L+M+1
VLM(:,LM,NIP) = VLMFULL(:,LM,NI)
enddo
enddo
enddo ion2
CALLMPI( M_sum_d(WDES%COMM_INTER, VLM, NDIM*LMMAX*WDES%NIONS) )
! test print
#ifdef MPI
if (GRIDC%COMM%NODE_ME == 1) then
#endif
open(unit=ofile, file='VLM.data', action='write')
do NR = 1, RNMAX
write(ofile, '(5E15.6)') PP%R%R(NR), VLM(NR,1,1), VLM(NR,2,1),VLM(NR,3,1), VLM(NR,4,1)
enddo
#ifdef MPI
endif
#endif
close(ofile)
#ifdef MPI
if (GRIDC%COMM%NODE_ME == 1) then
#endif
write(IU0,'(A)') ' External potential tranformed to VLM.'
#ifdef MPI
endif
#endif
deallocate(VLMFULL, CVEXT)
return
END SUBROUTINE EXTPT_CALC_VLM
!*************************** subroutine SETYLM_FINE_GRID **************************
! this subroutine setup the fine grid for the real space projection integeration
! I made up this layout to parallel the real space projection code
! Use LPAR to choose whether to distritute the local fine grid points or not
! Note:
! In CALC_VLM, we can not parallelize over uniform grid (GRIDC), so set LPAR = .Ture.
! to parallelize over the local fine grid
! In DE_DVRN, parallelize over GRIDC, so set LPAR = .False., making the local
! fine grid layout serial.
!**********************************************************************************
subroutine SETYLM_FINE_GRID( GRIDC, NGRID, YLM, LYMAX, LMMAX, &
XS, YS, ZS, DIST, IRDMAX, INDMAX, &
LPAR )
use prec
use constant
use mgrid
use mpimy
use asa
implicit none
type (grid_3d) GRIDC
integer :: NG, LYMAX, LMMAX, IRDMAX, NGRID
integer :: INDMAX
real(q) :: XS(IRDMAX), YS(IRDMAX), ZS(IRDMAX)
real(q) :: YLM(IRDMAX,LMMAX), DIST(IRDMAX)
integer :: ix, iy, iz, NODE_ME, NCPU, IND
logical :: LPAR
#ifdef MPI
NCPU = GRIDC%COMM%NCPU
NODE_ME = GRIDC%COMM%NODE_ME
#else
NCPU = 1
NODE_ME = 1
#endif
NG = 0
IND = 0
do ix = 1, NGRID
do iy = 1, NGRID
do iz = 1, NGRID
NG = NG + 1
if ( (.not. LPAR) .or. (mod(NG,NCPU)+1 == NODE_ME) ) then
IND = IND + 1
XS(IND) = (ix-0.5_q)*2.0_q/NGRID - 1.0_q
YS(IND) = (iy-0.5_q)*2.0_q/NGRID - 1.0_q
ZS(IND) = (iz-0.5_q)*2.0_q/NGRID - 1.0_q
DIST(IND) = sqrt(XS(IND)**2 + YS(IND)**2 + ZS(IND)**2)
if ( DIST(IND) > 1.0e-8 ) then
XS(IND) = XS(IND) / DIST(IND)
YS(IND) = YS(IND) / DIST(IND)
ZS(IND) = ZS(IND) / DIST(IND)
endif
endif
enddo
enddo
enddo
INDMAX = IND
call SETYLM(LYMAX,INDMAX,YLM,XS,YS,ZS)
return
end subroutine SETYLM_FINE_GRID
!*************************** function bspline **************************
! this function computes the interpolate value, using the interpolation
! coefficeint grid
! r should be the direct coordinate
!***********************************************************************
function bspline( r, s_grid, nx, ny, nz )
use prec
use constant
implicit none
integer :: nx, ny, nz
real(q) :: r(3), s_grid(nx, ny, nz), bspline
real(q) :: Bx, By, Bz, x, y, z, dx, dy, dz, h
integer :: dix, diy, diz, ix0, iy0, iz0, ix, iy, iz
bspline = 0.0_q
h = 1.0_q
! r ~ 0.0-1.0
ix0 = floor(r(1)*nx)+1
iy0 = floor(r(2)*ny)+1
iz0 = floor(r(3)*nz)+1
do diz = -1, 2
iz = mod(iz0+diz-1 + 10*nz, nz) + 1
dz = r(3)*nz - (iz0-1) - diz
Bz = Bfunc(dz, h)
do diy = -1, 2
iy = mod(iy0+diy-1 + 10*ny, ny) + 1
dy = r(2)*ny - (iy0-1) - diy
By = Bfunc(dy, h)
do dix = -1, 2
ix = mod(ix0+dix-1 + 10*nx, nx) + 1
dx = r(1)*nx - (ix0-1) - dix
Bx = Bfunc(dx, h)
bspline = bspline + s_grid(ix,iy,iz)*Bx*By*Bz
enddo
enddo
enddo
return
end function bspline
!*********************** SUBROUTINE BSPLINE_INTERPOLATE_COEFF ******************
! This subroutine conducts a 3-d B-spline interpolation
! f_grid: the actual input grid
! s_grid: the output spline coefficient grid
! sx, sy, sz : the 1-d coefficients for Delta-grid
!*******************************************************************************
subroutine BSPLINE_INTERPOLATE_COEFF( f_grid, s_grid, nx, ny, nz, di_smear, &
sx, sy, sz, GRIDC )
use prec
use constant
use mgrid
use mpimy
implicit none
integer :: nx, ny, nz, di_smear, ic, ir
real(q) :: f_grid(nx,ny,nz), s_grid(nx,ny,nz)
integer :: ix, iy, iz, ixp, iyp, izp, dix, diy, diz
real(q) :: sx(-di_smear:di_smear),sy(-di_smear:di_smear),sz(-di_smear:di_smear)
type (grid_3d) GRIDC
s_grid = 0.0_q
do ic = 1, GRIDC%RL%NCOL
if (GRIDC%RL%NFAST==3) then ! mpi version: x-> N2, y-> N3, z-> N1
ix = GRIDC%RL%I2(ic)
iy = GRIDC%RL%I3(ic)
else ! conventional version: x-> N1, y-> N2, z-> N3
iy = GRIDC%RL%I2(ic)
iz = GRIDC%RL%I3(ic)
endif
do ir = 1, GRIDC%RL%NROW
if (GRIDC%RL%NFAST==3) then
iz = ir
else
ix = ir
endif
do diz = -di_smear, di_smear
izp = mod(iz+diz + 10*nz-1, nz) + 1
do diy = -di_smear, di_smear
iyp = mod(iy+diy + 10*ny-1, ny) + 1
do dix = -di_smear, di_smear
ixp = mod(ix+dix + 10*nx-1, nx) + 1
s_grid(ixp,iyp,izp) = s_grid(ixp,iyp,izp) + &
sx(dix)*sy(diy)*sz(diz)*f_grid(ix,iy,iz)
enddo
enddo
enddo
enddo
enddo
CALLMPI( M_sum_d(GRIDC%COMM,s_grid,nx*ny*nz))
return
end subroutine BSPLINE_INTERPOLATE_COEFF
!*********************** SUBROUTINE CALC_Sdelta_COEFF ******************
! Calculate the 1-D B-spline coefficients for a "delta" grid
!***********************************************************************
subroutine CALC_Sdelta_COEFF(N, s, di_smear)
use prec
use constant
implicit none
integer :: N, di_smear
real(q) :: s(-di_smear:di_smear)
real(q) :: A(N,N), B(N)
integer :: i, j, ip, info
A = 0.0_q
B = 0.0_q
do i = 1, N
do j = -1, 1
ip = i+j
ip = mod( (ip-1 + 10*N), N ) + 1
if ( j==0 ) then
A(i,ip) = 4.0_q
else
A(i,ip) = 1.0_q
endif
enddo
enddo
B(1) = 1.0_q
! use mkl function to find solution
call DPOTRF( 'U', N, A, N, info ) ! Cholesky decomposition
call DPOTRS( 'U', N, 1, A, N, B, N, info ) ! solve Ax=B
do i = 1, di_smear
s(i) = B(i+1)
s(-i) = B(n-i+1)
enddo
s(0) = B(1)
s = s*6.0_q
return
end subroutine CALC_Sdelta_COEFF
!*********************** SUBROUTINE Bfunc ******************************
! This is the basis functions used in the B-spline interpolation
!***********************************************************************
function Bfunc(x, h)
use prec
use constant
implicit none
real(q) :: x, h, Bfunc
if ( x<-2*h .or. x>2*h) then
Bfunc = 0.0_q
else if ( x<-h .or. x>h) then
Bfunc = 1._q/6*(2*h-abs(x))**3
else
Bfunc = 2._q/3*h**3 - 0.5_q*x**2*(2*h-abs(x))
endif
return
end function Bfunc
!*********************** SUBROUTINE READ_ELEM **************************
! Wrappers to read elements from CVTOT
!***********************************************************************
FUNCTION READ_ELEM(GRIDC, CVTOT, I)
USE prec
USE base
USE lattice
USE mpimy
USE mgrid
USE poscar
USE constant
implicit none
TYPE (grid_3d) GRIDC
COMPLEX(q) CVTOT(DIMREAL(GRIDC%MPLWV))
INTEGER :: I
REAL(q) :: READ_ELEM
CALL READ_ELEM_(GRIDC, CVTOT, I, READ_ELEM)
return
END FUNCTION READ_ELEM
!*************************** FUNCTION RK *******************************
!
! This function calculates the radial projection functions
! Right now use the 1/r*sqrt(2/RCUT)*sin(pi*K/RCUT*r)
!
!***********************************************************************
FUNCTION RK( R, K, L )
use prec
use constant
implicit none
real(q) :: R, L, RK
integer :: K
if ( R < 0.0_q .or. R > L ) then
write(6,*) 'Internal Error: wrong R value for RK,', R, L
endif
if ( R > 1e-5 ) then
RK = sqrt(2._q/L)*sin(PI*K*R/L) / R
else
RK = sqrt(2._q/L)*PI*K/L
endif
return
END FUNCTION RK
!************************ FUNCTION FSWITCH *****************************
!
! This function calculates the switch function f(r) for V(r)
! In principle, can be anything as long as:
! f(RMAX) = 1
! f'(RMAX) = 0
! f(RMAX+Delta) = 0
!
!***********************************************************************
FUNCTION FSWITCH( R, RMAX, Delta )
use prec
use constant
implicit none
real*8 :: FSWITCH, R, RMAX, Delta
real*8 :: dr
if ( R<RMAX ) then
FSWITCH = 1.0
else if ( R>RMAX+Delta ) then
FSWITCH = 0.0
else
dr = (R - RMAX)/Delta
FSWITCH = cos(dr*PI/2.0)
endif
return
end function FSWITCH
!************************* ENERGY_DERIVATIVE ***************************
! This subroutine augments the electron density to give dE/dV(r)
! Theoretically it should be the AE density
! Pratically it should be softer than AECCAR2 but harder than CHGCAR
! Need:
! CRHODE: PAW occupancy matrix
! CHDEN: soft tilde(n) without compensation charge hat(n)
! defined in GRID_SOFT
! CHTOT: tilde(n) + hat(n) as input, correct dE/dV as output
!***********************************************************************
SUBROUTINE EXTPT_DE_DVRN(WDES, GRID_SOFT, GRIDC, SOFT_TO_C, LATT_CUR, &
P, T_INFO, LMDIM, CRHODE, CHTOT_, CHDEN, IRDMAX )
use prec
use base
use mgrid
use mpimy
use lattice
use poscar
use pseudo
use us
use paw
use wave
use asa
use constant
implicit none
TYPE (type_info) T_INFO
TYPE (potcar) P(T_INFO%NTYP)
TYPE (grid_3d) GRIDC
TYPE (grid_3d) GRID_SOFT
TYPE (latt) LATT_CUR
TYPE (transit) SOFT_TO_C
TYPE (wavedes) WDES
integer :: LMDIM, IRDMAX
COMPLEX(q) CHDEN(GRID_SOFT%MPLWV,WDES%NCDIJ)
OVERLAP CRHODE(LMDIM,LMDIM,WDES%NIONS,WDES%NCDIJ)
COMPLEX(q),TARGET :: CHTOT_(DIMREAL(GRIDC%MPLWV),WDES%NCDIJ)
! local variable
integer, external :: MAXL_AUG, MAXL1
integer :: ISP, NT, NI, NIP, LYMAX, LMMAX, INDMAX, IND, NDIM, I
integer :: LT, LPT, MT, MPT, IGRID
integer :: L, LP, LL, LLP, M, MP, MPLOW, FAKLLP, LM, LMP, RNMAX, FAKQ
integer :: LMINDX, ISTART, IEND, IC, NR, NC, NCPU, ICPU
integer :: NXMAX, NYMAX, NZMAX, NALLOC, NGF, JLL, JLM, ierror
integer, allocatable :: NLI(:), LLOCAL(:)
real(q), allocatable :: CWORK(:)
complex(q), allocatable :: CWORK2(:)
COMPLEX(q),POINTER :: CHTOT(:)
integer, parameter :: nprint=8
integer :: ofile
integer :: INDFMAX, INDF, K, ix, iy, iz, ixp, iyp, izp, dix, diy, diz
real(q) :: VLMK(KMAX), dV, RMAX, RCUT
real(q) :: RI(3), RVEC(3), RVECF(3), DRVEC(3), val
real(q) :: STMP1, STMP2, STMP3, STMP4, R, STMP1PREV, RPREV
real(q), allocatable :: DIST(:), QIJ(:)
real(q), allocatable :: XS(:),YS(:), ZS(:), SUM(:)
real(q), allocatable :: YLMF(:,:), XSF(:), YSF(:), ZSF(:), DISTF(:)
OVERLAP, allocatable :: CTMP(:,:,:,:)
character(100) :: fn
character(50) :: buffer
logical :: LINFO
! temparary storage matrix for rho_ij
ALLOCATE(CTMP(LMDIM,LMDIM,T_INFO%NIONS,WDES%NCDIJ))
allocate(CWORK(DIMREAL(GRIDC%MPLWV)), CWORK2((DIMREAL(GRIDC%MPLWV))))
LYMAX=MAXL_AUG(T_INFO%NTYP,P)
LMMAX=(LYMAX+1)**2 ! number of lm pairs
#ifdef MPI
NCPU = GRIDC%COMM%NCPU
#else
NCPU = 1
#endif
! Construct the full rho_ij matrix
CTMP=0
DO ISP=1,WDES%NCDIJ
DO NI=1,T_INFO%NIONS
NIP=NI_LOCAL(NI, WDES%COMM_INB)
IF (NIP/=0) THEN
CTMP(:,:,NI,ISP)=CRHODE(:,:,NIP,ISP)
ENDIF
ENDDO
ENDDO
! Reduce rho_ij for all nodes, assume realmode
CALLMPI( M_sum_d(WDES%COMM_INB,CTMP,LMDIM*LMDIM*T_INFO%NIONS*WDES%NCDIJ))
ALLOCATE( XS(IRDMAX), YS(IRDMAX), ZS(IRDMAX), &
DIST(IRDMAX),SUM(IRDMAX), &
NLI(IRDMAX), LLOCAL(IRDMAX))
NDIM=0
DO NT=1,T_INFO%NTYP
IF (ASSOCIATED(P(NT)%QPAW)) THEN
NDIM=MAX(NDIM, P(NT)%R%NMAX)
END IF
ENDDO
allocate(QIJ(NDIM))
#ifdef MPI
ICPU = WDES%COMM%NODE_ME
#else
ICPU = 1
#endif
NXMAX = EXTPT%nx
NYMAX = EXTPT%ny
NZMAX = EXTPT%nz
NALLOC = NGRID*NGRID*NGRID
allocate( YLMF(NALLOC,LMMAX), DISTF(NALLOC), XSF(NALLOC), &
& YSF(NALLOC), ZSF(NALLOC) )
CALLMPI( MPI_barrier( GRIDC%COMM, ierror ))
inquire(DIRECTORY="DERINFO", EXIST=LINFO) ! does external information exists?
CALLMPI( MPI_barrier( GRIDC%COMM, ierror ))
if ( .not. LINFO ) call system('mkdir DERINFO')
CALLMPI( MPI_barrier( GRIDC%COMM, ierror ))
write( buffer, "(I)" ) ICPU
buffer = trim(adjustl(buffer))
write( fn, "('DERINFO/',A)" ) buffer
ofile = 666+ICPU-1
! if DERINFO exists, open files to read
if ( LINFO ) then
if ( ICPU == 1 ) print *, 'read DERINFO'
open( unit=ofile, file=fn, action='read' )
! else, open files to write and prepare for the real calculations
else
if ( ICPU == 1 ) print *, 'write DERINFO'
open( unit=ofile, file=fn, action='write' )
! set up YLM for the fine integration grid, this is shared by all atom
! centers
call SETYLM_FINE_GRID( GRIDC, NGRID, YLMF, LYMAX, LMMAX, XSF &
, YSF, ZSF, DISTF, NALLOC, INDFMAX, .false.)
endif
! Do the actual work, only work with the total (alpha+beta) component
spin: do ISP = 1, 1 !WDES%NCDIJ
CHTOT => CHTOT_(:,ISP)
ion: do NI = 1, T_INFO%NIONS
if ( ICPU == 1 ) write(*,*) 'Ion:', NI
CWORK = 0.0_q
NT = T_INFO%ITYP(NI)
RNMAX = P(NT)%R%NMAX
RMAX = P(NT)%R%R(RNMAX)
RCUT = RMAX + Delta
dV = (RCUT*2/NGRID)**3
RI = T_INFO%POSION(1:3,NI)
! Find the relevant points on GRIDC, in my equation, there is no
! YLM(hat(r)), hence no need to use SETYLM_AUG
! In this version, we parallelize over local fine grid, so for each node,
! should have a complete grid point list, rather than local list
CALL SETRN(GRIDC,LATT_CUR,RI,P(NT)%PSDMAX*EXTPT%FCORR,NPSRNL,&
IRDMAX,INDMAX, 0.0,0.0,0.0,DIST(1),NLI(1),XS,YS,ZS,LLOCAL)
rn_loop: do IND = 1, INDMAX
SUM(IND) = 0.0_q
#ifdef MPI
if ( mod(IND,GRIDC%COMM%NCPU)+1 /= ICPU ) goto 321 ! if not local grid point, skip
#endif
SGRID = 0.0_q
if ( .not. LINFO ) then
RVEC(1) = XS(IND)*DIST(IND) ! in angstrom
RVEC(2) = YS(IND)*DIST(IND)
RVEC(3) = ZS(IND)*DIST(IND)
RVEC = matmul( transpose(LATT_CUR%B(1:3,1:3)), RVEC(1:3)) ! in direct
ix = nint((RVEC(1)+RI(1))*NXMAX)
iy = nint((RVEC(2)+RI(2))*NYMAX)
iz = nint((RVEC(3)+RI(3))*NZMAX)
ix = mod(ix + 10*NXMAX, NXMAX) + 1
iy = mod(iy + 10*NYMAX, NYMAX) + 1
iz = mod(iz + 10*NZMAX, NZMAX) + 1
! construct SGRID for a delta grid centered at (XS,YS,ZS)
do diz = -di_smear, di_smear
izp = mod(iz+diz + 10*NZMAX-1, NZMAX) + 1
do diy = -di_smear, di_smear
iyp = mod(iy+diy + 10*NYMAX-1, NYMAX) + 1
do dix = -di_smear, di_smear
ixp = mod(ix+dix + 10*NXMAX-1, NXMAX) + 1
SGRID(ixp,iyp,izp) = SGRID(ixp,iyp,izp) + &
SCX(dix)*SCY(diy)*SCZ(diz)
enddo
enddo
enddo
endif
! four DO loops over all the channels: l, m, lp, mp
LM = 1
l_loop: DO L =1,P(NT)%LMAX
LMP=LM
lp_loop: DO LP=L,P(NT)%LMAX
if ( .not. LINFO ) then
! setup QIJ
call SETQIJ(GRIDC, L,LP,P(NT),QIJ, NDIM)
endif
! Electron density scaled by occupation number
LL =P(NT)%LPS(L )
LLP=P(NT)%LPS(LP)
if ( .not. LINFO ) then
! setup the YLM3 lookup index
CALL YLM3LOOKUP(LL,LLP,LMINDX)
endif
m_loop: DO M=1,2*LL+1
MPLOW=1
IF (L==LP) MPLOW=M
mp_loop: DO MP=1,2*LLP+1 !MPLOW,2*LLP+1
FAKLLP=1
!IF (LMP+MP/=LM+M) FAKLLP=2
IF (L/=LP) FAKLLP=2
if ( .not. LINFO ) then
! setup YLM3 for each (l, m, lp, mp)
LMINDX=LMINDX+1
ISTART=INDCG(LMINDX)
IEND =INDCG(LMINDX+1)
endif
! Loop over L, M
! JL(IC) would point to L
! JS(IC) would point to L*(L+1) + M + 1, for YLM3(IC)
STMP3 = 0.0_q
if (LINFO) then
!*****************************
! read STMP3 from DERINFO
!*****************************
read(ofile,'(I4,I12,4I4,E25.10)') NIP,IGRID,LT,MT,LPT,MPT,STMP3
! check consistency
if (NIP/=NI .or. IGRID/=ABS(LLOCAL(IND)) .or. LT/=L .or. MT/=M &
.or. LPT/=LP .or. MPT/=MP ) then
write(6,*) 'Internal error, mismatch with DERINFO'
endif
else
!*****************************
! compute and write DERINFO
!*****************************
DO IC=ISTART,IEND-1
JLL = JL(IC)
JLM = JS(IC)
! find out VLMK for Delta(XS, YS, ZS)
VLMK(:) = 0.0_q
do K = 1, KMAX
! \sum YLM *RK *V(r) *dV, in fine grid
do INDF = 1, INDFMAX
R = DISTF(INDF)*RCUT ! in Angstrom
if ( R <= RCUT ) then
RVECF(1) = XSF(INDF)*R
RVECF(2) = YSF(INDF)*R
RVECF(3) = ZSF(INDF)*R
RVECF = matmul( transpose(LATT_CUR%B(1:3,1:3)), RVECF(1:3) )
DRVEC = RVECF - RVEC ! in direct
if(abs(DRVEC(1)*NXMAX)<di_smear+1 .and. &
abs(DRVEC(2)*NYMAX)<di_smear+1 .and. &
abs(DRVEC(3)*NZMAX)<di_smear+1) then
RVECF(:) = RVECF(:) + RI(:) ! grid coordinate in direct
val = bspline( RVECF, SGRID, NXMAX, NYMAX, NZMAX )
VLMK(K) = VLMK(K) + dV* FSWITCH(R,RMAX,Delta)* &
& YLMF(INDF,JLM)*RK(R,K,RCUT)*val
endif
endif
enddo
enddo
! loop over r
STMP2 = 0.0_q
do NR = RNMAX, 1, -1
R = P(NT)%R%R(NR)
if ( R>RMAX/KMAX .or. JLL==0 ) then
STMP1 = 0.0_q
do K = 1, KMAX
STMP1 = STMP1 + RK(R,K,RCUT) * VLMK(K)
enddo
STMP1PREV = STMP1
RPREV = R
else ! in the r->0 region
STMP1 = STMP1PREV * ((R/RPREV)**JLL)
endif
! \sum_r W(r)*QIJ(r)
STMP2 = STMP2 + STMP1*P(NT)%R%SI(NR)* &
(QIJ(NR) - P(NT)%QPAW(L,LP,JLL)*P(NT)%AUG(NR,JLL))
enddo
! \sum_LM C(LM, l, m, lp, mp)
STMP3 = STMP3 + STMP2*YLM3(IC)
ENDDO
! write the external info
write(ofile,'(I4,I12,4I4,E25.10)') NI,abs(LLOCAL(IND)),L,M,LP,MP,STMP3
endif
! \sum_llpmmp rho_ij
SUM(IND) = SUM(IND) + STMP3*CTMP(LM+M-1,LMP+MP-1,NI,ISP)*FAKLLP
ENDDO mp_loop
ENDDO m_loop
LMP=LMP+2*LLP+1
ENDDO lp_loop
LM =LM +2*LL +1
ENDDO l_loop
321 continue
enddo rn_loop
! synchronize the rn_loop
CALLMPI( M_sum_d(GRIDC%COMM, SUM, INDMAX) )
do IND = 1, INDMAX
! when adding, only deal with local grid points
if (LLOCAL(IND) > 0) then
CWORK(NLI(IND)) = CWORK(NLI(IND)) + SUM(IND)*NXMAX*NYMAX*NZMAX
endif
enddo
! FFT to reciprocal space
CALL FFT_RC_SCALE(CWORK(1),CWORK2(1),GRIDC)
CHTOT = CWORK2 + CHTOT
enddo ion
enddo spin
deallocate(XS,YS,ZS,DIST,SUM,NLI,LLOCAL,QIJ,YLMF,DISTF,XSF,YSF,ZSF,CTMP,CWORK)
close( ofile )
return
END SUBROUTINE EXTPT_DE_DVRN
end module mextpot
!************************ SUBROUTINE SETRN ****************************
!
! this subroutine performes the following tasks
! ) finds the points, which are within a certain cutoff around one ion
! ) calculates the distance of each of this points from the ion
! DISX,Y,Z are additional displacements of the ions
!
! mind that the cutoff-sphere extends up to PSDMAX*(NPSRNL-1)/NPSRNL
! This code is a direct copy from SETYLM_AUG, with the YLM part removed
!
!***********************************************************************
SUBROUTINE SETRN(GRID,LATT_CUR,POSION,PSDMAX,NPSRNL, &
IRMAX,INDMAX,DISX,DISY,DISZ,DIST,NLI, &
XS,YS,ZS, LLOCAL)
USE prec
USE mpimy
USE mgrid
USE lattice
USE asa
USE constant
IMPLICIT COMPLEX(q) (C)
IMPLICIT REAL(q) (A-B,D-H,O-Z)
TYPE (grid_3d) GRID
TYPE (latt) LATT_CUR
DIMENSION POSION(3)
DIMENSION DIST(IRMAX)
! work-arrays
DIMENSION NLI(IRMAX)
REAL(q) XS(IRMAX),YS(IRMAX),ZS(IRMAX)
INTEGER :: LLOCAL(IRMAX)
XS=0;YS=0;ZS=0; LLOCAL=0
!=======================================================================
! find lattice points contained within the cutoff-sphere
! mind that the cutoff-sphere extends up to PSDMAX*(NPSRNL-1)/NPSRNL
! which is a somewhat strange convention
!=======================================================================
F1=1._q/GRID%NGX
F2=1._q/GRID%NGY
F3=1._q/GRID%NGZ
ARGSC=NPSRNL/PSDMAX
!-----------------------------------------------------------------------
! restrict loop to points contained within a cubus around the ion
!-----------------------------------------------------------------------
D1= PSDMAX*LATT_CUR%BNORM(1)*GRID%NGX
D2= PSDMAX*LATT_CUR%BNORM(2)*GRID%NGY
D3= PSDMAX*LATT_CUR%BNORM(3)*GRID%NGZ
N3LOW= INT(POSION(3)*GRID%NGZ-D3+10*GRID%NGZ+.99_q)-10*GRID%NGZ
N2LOW= INT(POSION(2)*GRID%NGY-D2+10*GRID%NGY+.99_q)-10*GRID%NGY
N1LOW= INT(POSION(1)*GRID%NGX-D1+10*GRID%NGX+.99_q)-10*GRID%NGX
N3HI = INT(POSION(3)*GRID%NGZ+D3+10*GRID%NGZ)-10*GRID%NGZ
N2HI = INT(POSION(2)*GRID%NGY+D2+10*GRID%NGY)-10*GRID%NGY
N1HI = INT(POSION(1)*GRID%NGX+D1+10*GRID%NGX)-10*GRID%NGX
!-----------------------------------------------------------------------