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