1! --- 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt . 6! See Docs/Contributors.txt for a list of contributors. 7! --- 8!!@LICENSE 9 10 MODULE m_chkgmx 11 12! Used module procedures 13 use sys, only: die ! Termination routine 14 use m_minvec, only: minvec ! Finds a lattice basis of minimal length 15 use cellsubs, only: reclat ! Finds reciprocal unit cell 16 17! Used module parameters and variables 18 use precision, only: dp ! Double precision real kind 19! use parallel, only: Node ! My processor index 20 21 implicit none 22 23! Public procedures provided: 24 PUBLIC:: 25 . chkgmx, ! Checks that a given cutoff is consistent with a mesh 26 . meshKcut ! Returns the planewave cutoff of a periodic mesh 27 28! Public parameters, variables, and arrays: 29! none 30 31 PRIVATE ! Nothing is declared public beyond this point 32 33 CONTAINS 34!--------------------------------------------------------------------------- 35 36 SUBROUTINE CHKGMX (K,BG,MESH,G2MAX) 37 38 implicit none 39 real(dp),intent(in) :: K(3) ! Vector in reciprocal unit cell (BZ) 40 real(dp),intent(in) :: BG(3,3) ! Reciprocal lattice vectors BG(ixyz,ivec) 41 integer, intent(in) :: MESH(3) ! Number of mesh point divisions of each 42 ! real-space lattice vector. 43 real(dp),intent(inout):: G2MAX ! Planewave cutoff, which is reduced to 44 ! the maximum value supported by the mesh 45 ! (without aliasing), if this is smaller. 46 47 real(dp), parameter :: ZERO=0.0_dp,HALF=0.5_dp, 48 . TINY=1.0e-8_dp,BIG=1.0e20_dp 49 50 integer :: i, j, i1, i2, i3 51 real(dp):: bm(3,3), baux(3,3), ctransf(3,3), 52 . g(3), gmod, gmax, g2msh, r 53 54 DO I=1,3 55 DO J=1,3 56 BM(J,I)=BG(J,I)*MESH(I) 57 Baux(J,I)=BG(J,I)*MESH(I) 58 ENDDO 59 ENDDO 60! 61! Use Baux to avoid aliasing of arguments, as one is intent(in) 62! and the other intent(out)... 63! 64 CALL MINVEC (Baux,BM, ctransf) 65 GMAX=BIG 66 DO I3=-1,1 67 DO I2=-1,1 68 DO I1=-1,1 69 IF (I1.EQ.0.AND.I2.EQ.0.AND.I3.EQ.0) GO TO 20 70 G(1)=BM(1,1)*I1+BM(1,2)*I2+BM(1,3)*I3 71 G(2)=BM(2,1)*I1+BM(2,2)*I2+BM(2,3)*I3 72 G(3)=BM(3,1)*I1+BM(3,2)*I2+BM(3,3)*I3 73 GMOD=SQRT(SUM(G**2)) 74 R=HALF*GMOD-SUM(K*G)/GMOD 75 GMAX=MIN(GMAX,R) 76 20 CONTINUE 77 ENDDO 78 ENDDO 79 ENDDO 80 IF (GMAX.LT.ZERO) call die('CHKGMX: K NOT IN FIRST BZ') 81 G2MSH=GMAX*GMAX-TINY 82 IF (G2MSH.LT.G2MAX) THEN 83! if (Node.eq.0) then 84! WRITE(6,*) 'CHKGMX: MESH TOO SPARSE. GMAX HAS BEEN REDUCED' 85! WRITE(6,*) 'CHKGMX: OLD G2MAX =',G2MAX 86! WRITE(6,*) 'CHKGMX: NEW G2MAX =',G2MSH 87! endif 88 G2MAX=G2MSH 89 ENDIF 90 END SUBROUTINE CHKGMX 91 92!---------------------------------------------------------------------- 93 94 FUNCTION meshKcut( cell, nMesh ) 95 96 ! Finds mesh planewave cutoff 97 98 use precision, only: dp 99 100 implicit none 101 real(dp),intent(in):: cell(3,3) ! Unit cell vectors cell(iXYZ,iVector) 102 integer, intent(in):: nMesh(3) ! Mesh divisions of each vector 103 real(dp) :: meshKcut ! Mesh wavevector cutoff 104 105 real(dp):: g2max, k0(3), kcell(3,3) 106 107 ! Find reciprocal cell vectors 108 call reclat( cell, kcell, 1 ) 109 110 ! Initialize arguments to call chkgmx 111 k0(1:3) = 0._dp 112 g2max = huge(g2max) 113 114 ! chkgmx will reduce g2max to square of mesh cutoff 115 call chkgmx( k0, kcell, nMesh, g2max ) 116 117 ! Return wavevector cutoff 118 meshKcut = sqrt(g2max) 119 120 END FUNCTION meshKcut 121 122 END MODULE m_chkgmx 123 124