1!
2! Copyright (C) 2002-2011 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!------------------------------------------------------------------------------!
9  MODULE cell_nose
10!------------------------------------------------------------------------------!
11
12      USE kinds, ONLY : DP
13!
14      IMPLICIT NONE
15      SAVE
16
17      REAL(DP) :: xnhh0(3,3) = 0.0_DP
18      REAL(DP) :: xnhhm(3,3) = 0.0_DP
19      REAL(DP) :: xnhhp(3,3) = 0.0_DP
20      REAL(DP) :: vnhh(3,3)  = 0.0_DP
21      REAL(DP) :: temph      = 0.0_DP  !  Thermostat temperature (from input)
22      REAL(DP) :: fnoseh     = 0.0_DP  !  Thermostat frequency (from input)
23      REAL(DP) :: qnh        = 0.0_DP  !  Thermostat mass (computed)
24
25CONTAINS
26
27  subroutine cell_nose_init( temph_init, fnoseh_init )
28     USE constants, ONLY: pi, au_terahertz, k_boltzmann_au
29     REAL(DP), INTENT(IN) :: temph_init, fnoseh_init
30     ! set thermostat parameter for cell
31     qnh    = 0.0_DP
32     temph  = temph_init
33     fnoseh = fnoseh_init
34     if( fnoseh > 0.0_DP ) qnh = 2.0_DP * ( 3 * 3 ) * temph * k_boltzmann_au / &
35          (fnoseh*(2.0_DP*pi)*au_terahertz)**2
36    return
37  end subroutine cell_nose_init
38
39  subroutine cell_nosezero( vnhh, xnhh0, xnhhm )
40    real(DP), intent(out) :: vnhh(3,3), xnhh0(3,3), xnhhm(3,3)
41    xnhh0=0.0_DP
42    xnhhm=0.0_DP
43    vnhh =0.0_DP
44    return
45  end subroutine cell_nosezero
46
47  subroutine cell_nosevel( vnhh, xnhh0, xnhhm, delt )
48    implicit none
49    REAL(DP), intent(inout) :: vnhh(3,3)
50    REAL(DP), intent(in) :: xnhh0(3,3), xnhhm(3,3), delt
51    vnhh(:,:)=2.0_DP*(xnhh0(:,:)-xnhhm(:,:))/delt-vnhh(:,:)
52    return
53  end subroutine cell_nosevel
54
55  subroutine cell_noseupd( xnhhp, xnhh0, xnhhm, delt, qnh, temphh, temph, vnhh )
56    use constants, only: k_boltzmann_au
57    implicit none
58    REAL(DP), intent(out) :: xnhhp(3,3), vnhh(3,3)
59    REAL(DP), intent(in) :: xnhh0(3,3), xnhhm(3,3), delt, qnh, temphh(3,3), temph
60    integer :: i, j
61    do j=1,3
62      do i=1,3
63        xnhhp(i,j) = 2.0_DP*xnhh0(i,j)-xnhhm(i,j) + &
64             (delt**2/qnh)* k_boltzmann_au * (temphh(i,j)-temph)
65        vnhh(i,j) =(xnhhp(i,j)-xnhhm(i,j))/( 2.0_DP * delt )
66      end do
67    end do
68    return
69  end subroutine cell_noseupd
70
71
72  REAL(DP) function cell_nose_nrg( qnh, xnhh0, vnhh, temph, iforceh )
73    use constants, only: k_boltzmann_au
74    implicit none
75    REAL(DP) :: qnh, vnhh( 3, 3 ), temph, xnhh0( 3, 3 )
76    integer :: iforceh( 3, 3 )
77    integer :: i, j
78    REAL(DP) :: enij
79    cell_nose_nrg = 0.0_DP
80    do i=1,3
81      do j=1,3
82        enij = 0.5_DP*qnh*vnhh(i,j)*vnhh(i,j)+temph*k_boltzmann_au*xnhh0(i,j)
83        cell_nose_nrg = cell_nose_nrg + iforceh( i, j ) * enij
84      enddo
85    enddo
86    return
87  end function cell_nose_nrg
88
89  subroutine cell_nose_shiftvar( xnhhp, xnhh0, xnhhm )
90    !  shift values of nose variables to start a new step
91    implicit none
92    REAL(DP), intent(out) :: xnhhm(3,3)
93    REAL(DP), intent(inout) :: xnhh0(3,3)
94    REAL(DP), intent(in) :: xnhhp(3,3)
95      xnhhm = xnhh0
96      xnhh0 = xnhhp
97    return
98  end subroutine cell_nose_shiftvar
99
100
101  SUBROUTINE cell_nose_info ( delt )
102
103      use constants,     only: au_terahertz, pi
104      USE io_global,     ONLY: stdout
105      USE control_flags, ONLY: tnoseh
106
107      IMPLICIT NONE
108
109      REAL(DP), INTENT (IN) :: delt
110
111      INTEGER   :: nsvar
112      REAL(DP) :: wnoseh
113
114      IF( tnoseh ) THEN
115        !
116        IF( fnoseh <= 0.0_DP) &
117          CALL errore(' cell_nose_info ', ' fnoseh less than zero ', 1)
118        IF( delt <= 0.0_DP) &
119          CALL errore(' cell_nose_info ', ' delt less than zero ', 1)
120
121        wnoseh = fnoseh * ( 2.0_DP * pi ) * au_terahertz
122        nsvar  = ( 2.0_DP * pi ) / ( wnoseh * delt )
123
124        WRITE( stdout,563) temph, nsvar, fnoseh, qnh
125      END IF
126
127 563  format( //, &
128     &       3X,'cell dynamics with nose` temperature control:', /, &
129     &       3X,'Kinetic energy required   = ', f10.5, ' (Kelvin) ', /, &
130     &       3X,'time steps per nose osc.  = ', i5, /, &
131     &       3X,'nose` frequency           = ', f10.3, ' (THz) ', /, &
132     &       3X,'nose` mass(es)            = ', 20(1X,f10.3),//)
133
134    RETURN
135  END SUBROUTINE cell_nose_info
136
137
138!
139!------------------------------------------------------------------------------!
140   END MODULE cell_nose
141!------------------------------------------------------------------------------!
142
143