comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 21:55:16 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

kpsym

Browse source
This commit is contained in:
Frederick Stein 2026-07-17 14:54:44 +02:00
parent 3caaa28e3b
commit f52b992d46
1 changed files with 137 additions and 188 deletions
Download patch file Download diff file Expand all files Collapse all files

325
src/kpsym.F
View file

@ -170,6 +170,7 @@ CONTAINS
INTEGER :: i, ib(48), ib0(48), ihc, ihc0, ihg, ihg0, indpg, indpg0, invadd, istrin, iswght, &
isy, isy0, itype, j, k, l, li, li0, lmax, n, nc, nc0, ntot, ntvec0
INTEGER, DIMENSION(49, 1) :: f00
LOGICAL :: is_located
REAL(KIND=dp) :: a01(3), a02(3), a03(3), b01(3), b02(3), b03(3), b1(3), b2(3), b3(3), &
dtotstr, origin(3), origin0(3), proj1, proj2, proj3, r(3, 3, 48), r0(3, 3, 48), totstr, &
tvec0(3, 1), volum, vv0(3)
@ -197,32 +198,36 @@ CONTAINS
itype = 0
DO i = 1, nat
! Assign an atomic type (for internal purposes)
is_located = .FALSE.
IF (i /= 1) THEN
DO j = 1, (i - 1)
IF (ty(j) == ty(i)) THEN
! Type located
GOTO 178
is_located = .TRUE.
EXIT
END IF
END DO
! New type
END IF
itype = itype + 1
IF (itype > nsp) THEN
IF (iout > 0) THEN
WRITE (iout, '(A,I4,")")') &
' KPSYM| NUMBER OF ATOMIC TYPES EXCEEDS DIMENSION (NSP=)', &
nsp
IF (.NOT. is_located) THEN
itype = itype + 1
IF (itype > nsp) THEN
IF (iout > 0) THEN
WRITE (iout, '(A,I4,")")') &
' KPSYM| NUMBER OF ATOMIC TYPES EXCEEDS DIMENSION (NSP=)', &
nsp
END IF
IF (iout > 0) THEN
WRITE (iout, '(" KPSYM| THE ARRAY TY IS:",/,9(1X,10I7,/))') &
(ty(j), j=1, nat)
END IF
CPABORT('K290: FATAL ERROR')
END IF
ELSE
IF (iout > 0) THEN
WRITE (iout, '(" KPSYM| THE ARRAY TY IS:",/,9(1X,10I7,/))') &
(ty(j), j=1, nat)
WRITE (iout, '(" KPSYM|",6X,I5,I6,3F10.5)') &
i, ty(i), (xkapa(j, i), j=1, 3)
END IF
CALL stopgm('K290', 'FATAL ERROR')
END IF
178 CONTINUE
IF (iout > 0) THEN
WRITE (iout, '(" KPSYM|",6X,I5,I6,3F10.5)') &
i, ty(i), (xkapa(j, i), j=1, 3)
END IF
END DO
! ==--------------------------------------------------------------==
@ -436,7 +441,7 @@ CONTAINS
WRITE (iout, '(A,/,1X,3F10.5)') &
' KPSYM| CONSTANT VECTOR SHIFT (MACDONALD) OF THIS MESH:', wvk0
END IF
IF (ABS(iq1) + ABS(iq2) + ABS(iq3) == 0) GOTO 710
IF (ABS(iq1) + ABS(iq2) + ABS(iq3) == 0) RETURN
IF (ABS(istriz) /= 1) THEN
IF (iout > 0) THEN
WRITE (iout, '(" KPSYM| INVALID SWITCH FOR SYMMETRIZATION",I10)') istriz
@ -444,7 +449,7 @@ CONTAINS
IF (iout > 0) THEN
WRITE (iout, '(" KPSYM| INVALID SWITCH FOR SYMMETRIZATION",I10)') istriz
END IF
CALL stopgm('K290', 'ISTRIZ WRONG ARGUMENT')
CPABORT('K290: ISTRIZ WRONG ARGUMENT')
END IF
IF (iout > 0) THEN
WRITE (iout, '(" KPSYM| SYMMETRIZATION SWITCH: ",I3)', advance="no") istriz
@ -480,8 +485,7 @@ CONTAINS
IF (iout > 0) THEN
WRITE (iout, *) ' KPSYM| NO DUPLICATION FOUND'
END IF
CALL stopgm('ERROR', &
'SOMETHING IS WRONG IN GROUP DETERMINATION')
CPABORT('ERROR: SOMETHING IS WRONG IN GROUP DETERMINATION')
END IF
nc = nc0
DO i = 1, nc0
@ -517,7 +521,7 @@ CONTAINS
WRITE (iout, '(/," KPSYM|",1X,I5," SPECIAL POINTS GENERATED")') ntot
END IF
IF (ntot == 0) THEN
GOTO 710
RETURN
ELSE IF (ntot < 0) THEN
IF (iout > 0) THEN
WRITE (iout, '(A,I5,/,A,/,A)') ' KPSYM| DIMENSION NKPOINT =', nkpoint, &
@ -571,10 +575,6 @@ CONTAINS
IF (iout > 0) THEN
WRITE (iout, '(24X,"TOTAL:",I8)') iswght
END IF
! ==--------------------------------------------------------------==
710 CONTINUE
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE k290s
! **************************************************************************************************
@ -824,8 +824,6 @@ CONTAINS
(A(IL, JL)*A(IU, JU) - A(IL, JU)*A(IU, JL))
END DO
END DO
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE calbrec
! ==================================================================
! **************************************************************************************************
@ -895,7 +893,7 @@ CONTAINS
f0(49, i) = i
END DO
DO k2 = 2, nat
IF (ty(1) /= ty(k2)) GOTO 100
IF (ty(1) /= ty(k2)) CYCLE
DO i = 1, 3
xb(i) = x(i, k2) - x(i, 1)
END DO
@ -914,7 +912,6 @@ CONTAINS
tvec(i, ntvec) = vr(i)
END DO
END IF
100 CONTINUE
END DO
! ==-------------------------------------------------------------==
DO i = 1, 3
@ -950,15 +947,12 @@ CONTAINS
END DO
! Calculate new API
CALL calbrec(ap, api)
GOTO 200 ! EXIT
EXIT
END IF
END IF
END DO
200 CONTINUE
END DO
END IF
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE primlatt
! ==================================================================
! **************************************************************************************************
@ -1026,7 +1020,7 @@ CONTAINS
nc = 0
! Constructs rotation operations.
CALL rot1(ihc, r)
DO n = 1, nr
rotation_operations: DO n = 1, nr
ib(n) = 0
! Rotate the A1,2,3 vectors by rotation No. N
DO k = 1, 3
@ -1042,12 +1036,11 @@ CONTAINS
tr = tr + ABS(vr(i))
END DO
! If VR.ne.0, then XA cannot be a multiple of a lattice vector
IF (tr > delta) GOTO 140
IF (tr > delta) CYCLE rotation_operations
END DO
nc = nc + 1
ib(nc) = n
140 CONTINUE
END DO
END DO rotation_operations
! ==------------------------------------------------------------==
! IHG stands for holohedral group number.
IF (ihc == 0) THEN
@ -1066,8 +1059,6 @@ CONTAINS
RETURN
END IF
END DO
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE pgl1
! ==================================================================
! **************************************************************************************************
@ -1120,8 +1111,6 @@ CONTAINS
! VR(I) = - MOD(real(VR(I),kind=dp),1._dp)
vr(i) = NINT(vr(i)) - vr(i)
END DO
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE rlv3
! ==================================================================
! **************************************************************************************************
@ -1158,7 +1147,7 @@ CONTAINS
! == THE ATOM TRANSFORMATION TABLE,F0, ==
! == THE FRACTIONAL TRANSLATIONS,V, ==
! == ASSOCIATED WITH EACH ROTATION. ==
! == SUBROUTINES NEEDED: RLV3 CHECKRLV3 SYMMORPHIC STOPGM XSTRING ==
! == SUBROUTINES NEEDED: RLV3 CHECKRLV3 SYMMORPHIC XSTRING ==
! == MAY 14TH,1998: A LOT OF CHANGES (ARGUMENTS) ==
! == BETTER DETERMINATION OF V ==
! == SEP 15TH,1998: DETERMINATION OF FRACTIONAL TRANSLATIONAL VEC.==
@ -1276,36 +1265,31 @@ CONTAINS
! First we determine for VR=(/0,0,0/)
! IMPORTANT IF NOT UNIQUE ATOMS FOR DETERMINATION OF SYMMORPHIC
CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
IF (oksym) THEN
GOTO 190
END IF
! Now we try other possible VR
! F0(49,1:NAT) has only inequivalent atom indexes for translation
DO k2 = 1, nat
IF (f0(49, k2) < k2) GOTO 185
IF (ty(1) /= ty(k2)) GOTO 185
DO i = 1, 3
xb(i) = rx(i, 1) - x(i, k2)
IF (.NOT. oksym) THEN
! Now we try other possible VR
! F0(49,1:NAT) has only inequivalent atom indexes for translation
DO k2 = 1, nat
IF (f0(49, k2) < k2) CYCLE
IF (ty(1) /= ty(k2)) CYCLE
DO i = 1, 3
xb(i) = rx(i, 1) - x(i, k2)
END DO
! A translation vector VR is defined.
CALL rlv3(ai, xb, vr, il, delta)
! ==----------------------------------------------------------==
! == SUBROUTINE RLV3 REMOVES A DIRECT LATTICE VECTOR FROM XB ==
! == LEAVING THE REMAINDER IN VR. IF A NONZERO LATTICE ==
! == VECTOR WAS REMOVED, IL IS MADE NONZERO. ==
! == VR STANDS FOR V-REFERENCE. ==
! == VR IS NOT GIVEN IN CARTESIAN COORDINATES BUT ==
! == IN THE SYSTEM A1,A2,A3. K.K., 23.10.1979 ==
! ==----------------------------------------------------------==
CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
IF (oksym) EXIT
END DO
! A translation vector VR is defined.
CALL rlv3(ai, xb, vr, il, delta)
! ==----------------------------------------------------------==
! == SUBROUTINE RLV3 REMOVES A DIRECT LATTICE VECTOR FROM XB ==
! == LEAVING THE REMAINDER IN VR. IF A NONZERO LATTICE ==
! == VECTOR WAS REMOVED, IL IS MADE NONZERO. ==
! == VR STANDS FOR V-REFERENCE. ==
! == VR IS NOT GIVEN IN CARTESIAN COORDINATES BUT ==
! == IN THE SYSTEM A1,A2,A3. K.K., 23.10.1979 ==
! ==----------------------------------------------------------==
CALL checkrlv3(n, nat, ty, rx, x, vr, f0, ai, isc, nodupli, oksym, delta)
IF (oksym) THEN
GOTO 190
END IF
185 CONTINUE
END DO
iis(l) = 0
GOTO 210
190 CONTINUE
IF (.NOT. oksym) iis(l) = 0
CYCLE
END IF
nca = nca + 1
DO i = 1, 3
v(i, nca) = vr(i)
@ -1317,7 +1301,6 @@ CONTAINS
! == GIVEN IN THE SYSTEM A1,A2,A3. ==
! == K.K., 23.10. 1979 ==
! ==------------------------------------------------------------==
210 CONTINUE
END DO
! Remove unused operations
i = 0
@ -1326,14 +1309,13 @@ CONTAINS
li = 0
DO n = 1, nc
l = ib(n)
IF (iis(l) == 0) GOTO 230 ! CYCLE
IF (iis(l) == 0) CYCLE
i = i + 1
ib(i) = ib(n)
IF (ib(i) == ni) li = i
DO k = 1, nat
f0(i, k) = f0(n, k)
END DO
230 CONTINUE
END DO
! ==--------------------------------------------------------------==
nc = i
@ -1357,7 +1339,7 @@ CONTAINS
IF (iout > 0) THEN
WRITE (iout, '(" ATFTM1! IHG=",A," NC=",I2)') icst(ihg), nC
END IF
CALL stopgm('ATFTM1', 'NUMBER OF ROTATION NULL')
CPABORT('ATFTM1: NUMBER OF ROTATION NULL')
! Triclinic system
ELSE IF (nc == 1) THEN
! IB=1
@ -1501,7 +1483,7 @@ CONTAINS
IF (iout > 0) THEN
WRITE (iout, '(" ATFTM1! IHG=",A," NC=",I2)') icst(ihg), nC
END IF
CALL stopgm('ATFTM1', 'NUMBER OF ROTATION NULL')
CPABORT('ATFTM1: NUMBER OF ROTATION NULL')
! Triclinic system
ELSE IF (nc == 1) THEN
! IB=1
@ -1766,7 +1748,7 @@ CONTAINS
isc(ia) = 0
END DO
! Now we check if ROT(N)+VR gives a correct symmetry.
DO ia = 1, nat
atom: DO ia = 1, nat
DO ib = 1, nat
IF (ty(ia) == ty(ib) .AND. isc(ib) == 0) THEN
xb(1) = rx(1, ia) - x(1, ib)
@ -1782,16 +1764,13 @@ CONTAINS
f0(n, ia) = ib
! IR+VR is the good one: another symmetry operation
! Next atom
GOTO 100
CYCLE atom
END IF
END IF
END DO
! VR is not the correct translation vector
RETURN
100 CONTINUE
END DO
! ==--------------------------------------------------------------==
RETURN
END DO atom
END SUBROUTINE checkrlv3
! ==================================================================
! **************************************************************************************************
@ -1985,8 +1964,6 @@ CONTAINS
RETURN
END IF
END DO
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE symmorphic
! ==================================================================
! **************************************************************************************************
@ -2147,8 +2124,6 @@ CONTAINS
r(3, 3, nv) = -r(3, 3, n)
END DO
END IF
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE rot1
! ==================================================================
! **************************************************************************************************
@ -2388,7 +2363,18 @@ CONTAINS
END DO
END DO
! Check that WVA is inside the 1 Bz.
IF (.NOT. inside_bz(wva, rsdir, nplane, delta)) GOTO 450
IF (.NOT. inside_bz(wva, rsdir, nplane, delta)) THEN
IF (iout > 0) THEN
WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
END IF
IF (iout > 0) THEN
WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
' THE VECTOR ', wva, &
' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
' BY ROTATION NO. ', ibrav(iop), ' IS OUTSIDE THE 1BZ'
END IF
CPABORT('SPPT2: VECTOR OUTSIDE THE 1BZ')
END IF
! Place WVA in list
iplace = 0
CALL mesh(iout, wva, iplace, igarb0, igarbg, &
@ -2396,7 +2382,15 @@ CONTAINS
! If WVA was new (and therefore inserted),
! IPLACE is the number.
IF (iplace > 0) imesh = iplace
IF (iplace > nkpoint) GOTO 470
IF (iplace > nkpoint) THEN
IF (iout > 0) THEN
WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
END IF
IF (iout > 0) THEN
WRITE (iout, *) 'MESH SIZE EXCEEDS NKPOINT=', nkpoint
END IF
CPABORT('SPPT2: MESH SIZE EXCEEDED')
END IF
END DO
ELSE
! Place WVK in list
@ -2404,7 +2398,15 @@ CONTAINS
CALL mesh(iout, wvk, iplace, igarb0, igarbg, &
nkpoint, nhash, list, rlist, delta)
imesh = iplace
IF (iplace > nkpoint) GOTO 470
IF (iplace > nkpoint) THEN
IF (iout > 0) THEN
WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
END IF
IF (iout > 0) THEN
WRITE (iout, *) 'MESH SIZE EXCEEDS NKPOINT=', nkpoint
END IF
CPABORT('SPPT2: MESH SIZE EXCEEDED')
END IF
END IF
END DO
END DO
@ -2448,7 +2450,7 @@ CONTAINS
proja(3) = proja(3) + wva(k)*a3(k)
END DO
! Now loop over all the rest of the mesh points
DO j = (i + 1), imesh
mesh_point: DO j = (i + 1), imesh
jplace = j
CALL mesh(iout, wvk, jplace, igarb0, igarbg, &
nkpoint, nhash, list, rlist, delta)
@ -2464,15 +2466,14 @@ CONTAINS
! Check (PROJA - PROJB): Is it integral ?
DO k = 1, 3
diff = proja(k) - projb(k)
IF (ABS(REAL(NINT(diff), kind=dp) - diff) > delta) GOTO 280
IF (ABS(REAL(NINT(diff), kind=dp) - diff) > delta) CYCLE mesh_point
END DO
! DIFF is integral: remove WVK from mesh:
CALL remove(wvk, jplace, igarb0, igarbg, &
nkpoint, nhash, list, rlist, delta)
! If WVK actually removed, increment IREMOV
IF (jplace > 0) iremov = iremov + 1
280 CONTINUE
END DO
END DO mesh_point
END DO
IF (iremov > 0 .AND. iout > 0) THEN
WRITE (iout, '(A,A,/,A,1X,I6,A,/)') &
@ -2485,9 +2486,9 @@ CONTAINS
! == IN THE MESH OF WAVEVECTORS, NOW SEARCH FOR EQUIVALENT POINTS:==
! == THE INVERSION (TIME REVERSAL !) MAY BE USED. ==
! ==--------------------------------------------------------------==
DO iwvk = 1, imesh
! IF(INCLUD(IWVK) == YES) GOTO 350
IF (BTEST(includ(iwvk), 0)) GOTO 350
mesh_vector: DO iwvk = 1, imesh
! IF(INCLUD(IWVK) == YES) CYCLE mesh_vector
IF (BTEST(includ(iwvk), 0)) CYCLE mesh_vector
! IWVK has not been encountered previously: new special point,
! (only if WVK is not a garbage vector, however.)
! INCLUD(IWVK) = YES
@ -2498,7 +2499,7 @@ CONTAINS
! Find out whether Wvk is in the garbage list
CALL garbag(wvk, igarbage, igarb0, &
nkpoint, nhash, list, rlist, delta)
IF (igarbage > 0) GOTO 350
IF (igarbage > 0) CYCLE mesh_vector
ntot = ntot + 1
! Give the index in the special k points table.
includ(iwvk) = includ(iwvk) + ntot*2
@ -2508,7 +2509,7 @@ CONTAINS
lwght(ntot) = 1
! ==-----------------------------------------------------------==
! Find all the equivalent points (symmetry given by atoms)
DO n = 1, nc
group_element: DO n = 1, nc
! Rotate:
DO i = 1, 3
wva(i) = 0._dp
@ -2533,27 +2534,37 @@ CONTAINS
IF (igarbage == 0) THEN
! I think this case is impossible (NC <= NCBRAV)
! Error message
GOTO 490
IF (iout > 0) THEN
WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
END IF
IF (iout > 0) THEN
WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
' THE VECTOR ', wva, &
' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
' BY ROTATION NO. ', ib(n), ' IS NOT IN THE LIST'
END IF
CPABORT('SPPT2: VECTOR NOT IN THE LIST')
END IF
END IF
END IF
! Find out whether WVA is in the garbage list
CALL garbag(wva, igarbage, igarb0, &
nkpoint, nhash, list, rlist, delta)
IF (igarbage > 0) GOTO 370
IF (igarbage > 0) CYCLE group_element
! Was WVA encountered before ?
! IF(INCLUD(IPLACE) == YES) GOTO 364
IF (BTEST(includ(iplace), 0)) GOTO 364
! Increment weight.
lwght(ntot) = lwght(ntot) + 1
lrot(lwght(ntot), ntot) = ib(n)*ibsign
! INCLUD(IPLACE) = YES
includ(iplace) = IBSET(includ(iplace), 0)
! This k-point is an image of a special k-point.
! Put the index of the special k-point.
includ(iplace) = includ(iplace) + ntot*2
! IF(INCLUD(IPLACE) == YES) skip
IF (.NOT. BTEST(includ(iplace), 0)) THEN
! Increment weight.
lwght(ntot) = lwght(ntot) + 1
lrot(lwght(ntot), ntot) = ib(n)*ibsign
! INCLUD(IPLACE) = YES
includ(iplace) = IBSET(includ(iplace), 0)
! This k-point is an image of a special k-point.
! Put the index of the special k-point.
includ(iplace) = includ(iplace) + ntot*2
END IF
364 CONTINUE
IF (ibsign == -1 .OR. inv == 0) GOTO 370
IF (ibsign == -1 .OR. inv == 0) CYCLE group_element
! The case where we also apply the inversion to WVA
! Repeat the search, but for -WVA
ibsign = -1
@ -2561,10 +2572,8 @@ CONTAINS
wva(i) = -wva(i)
END DO
GOTO 363
370 CONTINUE
END DO
350 CONTINUE
END DO
END DO group_element
END DO mesh_vector
! ==--------------------------------------------------------------==
! == TOTAL NUMBER OF SPECIAL POINTS: NTOT ==
! == BEFORE USING THE LIST WVKL AS WAVE VECTORS, THEY HAVE TO BE ==
@ -2603,44 +2612,6 @@ CONTAINS
WRITE (iout, *)
END IF
END IF
RETURN
! ==--------------------------------------------------------------==
! == ERROR MESSAGES ==
! ==--------------------------------------------------------------==
450 CONTINUE
IF (iout > 0) THEN
WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
END IF
IF (iout > 0) THEN
WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
' THE VECTOR ', wva, &
' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
' BY ROTATION NO. ', ibrav(iop), ' IS OUTSIDE THE 1BZ'
END IF
CALL stopgm('SPPT2', 'VECTOR OUTSIDE THE 1BZ')
! ==--------------------------------------------------------------==
470 CONTINUE
IF (iout > 0) THEN
WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
END IF
IF (iout > 0) THEN
WRITE (iout, *) 'MESH SIZE EXCEEDS NKPOINT=', nkpoint
END IF
CALL stopgm('SPPT2', 'MESH SIZE EXCEEDED')
! ==--------------------------------------------------------------==
490 CONTINUE
IF (iout > 0) THEN
WRITE (iout, '(A,/)') ' SUBROUTINE SPPT2 *** FATAL ERROR ***'
END IF
IF (iout > 0) THEN
WRITE (iout, '(A,3F10.4,/,A,3F10.4,A,/,A,I3,A)') &
' THE VECTOR ', wva, &
' GENERATED FROM ', wvk, ' IN THE BASIC MESH', &
' BY ROTATION NO. ', ib(n), ' IS NOT IN THE LIST'
END IF
CALL stopgm('SPPT2', 'VECTOR NOT IN THE LIST')
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE sppt2
! **************************************************************************************************
!> \brief ...
@ -2729,7 +2700,7 @@ CONTAINS
' SUBROUTINE MESH *** FATAL ERROR *** LINKED LIST', &
' TOO LONG ***', ' CHOOSE A BETTER HASH-FUNCTION'
END IF
CALL stopgm('MESH', 'WARNING')
CPABORT('MESH: WARNING')
! WVK was not found
130 CONTINUE
IF (iplace == -1) THEN
@ -2748,7 +2719,7 @@ CONTAINS
' ISTORE=', istore, ' EXCEEDS DIMENSIONS', &
' WVK = ', wvk
END IF
CALL stopgm('MESH', 'WARNING')
CPABORT('MESH: WARNING')
END IF
list(istore) = nil
DO i = 1, 3
@ -2769,16 +2740,16 @@ CONTAINS
! == Return a wavevector (IPLACE > 0) ==
! ==--------------------------------------------------------------==
ipoint = iplace
IF (ipoint >= istore) GOTO 190
DO i = 1, 3
wvk(i) = rlist(i, ipoint)
END DO
RETURN
IF (ipoint < istore) THEN
DO i = 1, 3
wvk(i) = rlist(i, ipoint)
END DO
RETURN
END IF
END IF
! ==--------------------------------------------------------------==
! == Error - beyond list ==
! ==--------------------------------------------------------------==
190 CONTINUE
IF (iout > 0) THEN
WRITE (iout, '(A,/,A,I5,A,/)') &
' SUBROUTINE MESH *** WARNING ***', &
@ -2788,8 +2759,6 @@ CONTAINS
DO i = 1, 3
wvk(i) = 1.0e38_dp
END DO
! ==--------------------------------------------------------------==
RETURN
END SUBROUTINE mesh
! **************************************************************************************************
!> \brief ...
@ -2871,9 +2840,7 @@ CONTAINS
ipoint = list(ihash)
END DO
! List too long
CALL stopgm('MESH', 'LIST TOO LONG')
! ==--------------------------------------------------------------==
RETURN
CPABORT('MESH: LIST TOO LONG')
END SUBROUTINE remove
! **************************************************************************************************
!> \brief ...
@ -2937,9 +2904,7 @@ CONTAINS
ipoint = list(ihash)
END DO
! List too long
CALL stopgm('GARBAG', 'LIST TOO LONG')
! ==--------------------------------------------------------------==
RETURN
CPABORT('GARBAG: LIST TOO LONG')
END SUBROUTINE garbag
! **************************************************************************************************
@ -3101,9 +3066,9 @@ CONTAINS
i1 = nb1 + 1 - n1
DO n2 = 1, nnb2
i2 = nb2 + 1 - n2
DO n3 = 1, nnb3
last_dim: DO n3 = 1, nnb3
i3 = nb3 + 1 - n3
IF (i1 == 0 .AND. i2 == 0 .AND. i3 == 0) GOTO 150
IF (i1 == 0 .AND. i2 == 0 .AND. i3 == 0) CYCLE last_dim
DO i = 1, 3
bvec(i) = REAL(i1, kind=dp)*b1(i) + REAL(i2, kind=dp)*b2(i) + &
REAL(i3, kind=dp)*b3(i)
@ -3115,7 +3080,7 @@ CONTAINS
! 1/2*BVEC is outside the Bz - skip this direction
! The 1.e-6_dp takes care of single points touching the Bz,
! and of the -(plane)
IF (ABS(projct) > 0.5_dp - delta) GOTO 150
IF (ABS(projct) > 0.5_dp - delta) CYCLE last_dim
END DO
! 1/2*BVEC further confines the 1Bz - include into RSDIR
nplane = nplane + 1
@ -3125,8 +3090,7 @@ CONTAINS
END DO
! Length squared
rsdir(4, nplane) = bvec(1)**2 + bvec(2)**2 + bvec(3)**2
150 CONTINUE
END DO
END DO last_dim
END DO
END DO
@ -3140,19 +3104,4 @@ CONTAINS
END SUBROUTINE bzdefine
! **************************************************************************************************
!> \brief ...
!> \param a ...
!> \param b ...
! **************************************************************************************************
SUBROUTINE stopgm(a, b)
CHARACTER(LEN=*) :: a, b
CALL cp_warn(a, b)
CPABORT("stopgm@kpsym")
END SUBROUTINE stopgm
! **************************************************************************************************
END MODULE kpsym