From f52b992d46c7665bd2a8e8b7ee8a1d2b3a846e8d Mon Sep 17 00:00:00 2001 From: Frederick Stein Date: Fri, 17 Jul 2026 14:54:44 +0200 Subject: [PATCH] kpsym --- src/kpsym.F | 325 ++++++++++++++++++++++------------------------------ 1 file changed, 137 insertions(+), 188 deletions(-) diff --git a/src/kpsym.F b/src/kpsym.F index 4342c92f8b..57c87992b7 100644 --- a/src/kpsym.F +++ b/src/kpsym.F @@ -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