mirror of https://gitlab.com/QEF/q-e.git
Moving back to Modules directory prior to submission of Autopilot Suite
Here are former logs from CPV. 1.7 ballabio 2005-8-17 moved & in continued line to 6th position, ifort9 wants it [Gerardo] 1.6 kkudin 2005-07-29 Rescaled the target kinetic energy by ndega/(3*nat) for massive Nose chains, added a way to turn off the common thermostat on top of the massive Nose Kostya 1.5 kkudin 2005-07-29 For more than 1 Nose (chain) thermostat per system added a common thermostat on top of all the other ones Kostya 1.4 sbraccia 2005-07-18 Greneral cleanup. NEB works again also with the CP code. C.S 1.3 kkudin 2005-07-05 by Kostya This patch adds "massive" Nose-Hoover chains for ions (i.e. each ion can have a separate NH chain attached to it) Some fixes are still needed in different places: -the information on the number of NH chains [nhpdim] needs to be saved and read from the restart file (not done now) -the NH velocities also need to be all saved [nhpdim*nhpcl] -an input option needs to be added to zero out the NH velocities during a restart in order to permit "on the fly" thermostat changes -deallocation of the module variables gives glibc error with IFC 8.0 & 8.1 in cpr.f90 1.2 ballabio 2005-05-18 more end subroutine --> end subroutine name [Gerardo] 1.1 sbraccia 2005-05-16 ions_base splitted in three different files: ions_base.f90 (still in Modules), ions_positions.f90 and ions_nose (in CPV). C.S. git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@2114 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
parent
2680163daf
commit
c1f9da038e
|
@ -0,0 +1,376 @@
|
|||
!
|
||||
! Copyright (C) 2002-2005 FPMD-CPV groups
|
||||
! This file is distributed under the terms of the
|
||||
! GNU General Public License. See the file `License'
|
||||
! in the root directory of the present distribution,
|
||||
! or http://www.gnu.org/copyleft/gpl.txt .
|
||||
!
|
||||
!------------------------------------------------------------------------------!
|
||||
MODULE ions_nose
|
||||
!------------------------------------------------------------------------------!
|
||||
|
||||
USE kinds, ONLY: dbl
|
||||
USE parameters, ONLY: nhclm, natx
|
||||
!
|
||||
IMPLICIT NONE
|
||||
! Some comments are in order on how Nose-Hoover chains work here (K.N. Kudin)
|
||||
! the present code allows one to use "massive" Nose-Hoover chains:
|
||||
! TOBIAS DJ, MARTYNA GJ, KLEIN ML
|
||||
! JOURNAL OF PHYSICAL CHEMISTRY 97 (49): 12959-12966 DEC 9 1993
|
||||
!
|
||||
! currently "easy" input options allow one chain per atomic type (nhptyp=1)
|
||||
! and one chain per atom (nhptyp=2), but other options could be added too
|
||||
! one chain for the whole system is specified by nhptyp=0 (or nothing)
|
||||
!
|
||||
! nhpdim is the total number of the resulting NH chains
|
||||
! nhpend is 1 if there is a chain above all chains, otherwise it is 0
|
||||
! array atm2nhp(1:nat) gives the chain number from the atom list (which
|
||||
! is sorted by type)
|
||||
! anum2nhp is the number of degrees of freedom per chain (now just 3*nat_i)
|
||||
! ekin2nhp is the kinetic energy of the present chain
|
||||
! gkbt2nhp are the NH chain parameters
|
||||
! qnp are the chain masses, qnp_ is a temporary array for now
|
||||
! see subroutine ions_nose_allocate on what are the dimensions of these
|
||||
! variables
|
||||
! nhclm is now mostly not used, needs to be cleaned up at some point
|
||||
!
|
||||
INTEGER :: nhpcl, ndega, nhpdim, nhptyp, nhpend
|
||||
INTEGER :: atm2nhp(natx)
|
||||
INTEGER, ALLOCATABLE :: anum2nhp(:)
|
||||
REAL(dbl), ALLOCATABLE :: vnhp(:), xnhp0(:), xnhpm(:), xnhpp(:), &
|
||||
ekin2nhp(:), gkbt2nhp(:), qnp(:), qnp_(:)
|
||||
|
||||
REAL(dbl) :: gkbt = 0.0d0
|
||||
REAL(dbl) :: kbt = 0.0d0
|
||||
REAL(dbl) :: tempw = 0.0d0
|
||||
REAL(dbl) :: fnosep( nhclm ) = 0.0d0
|
||||
|
||||
!------------------------------------------------------------------------------!
|
||||
CONTAINS
|
||||
!------------------------------------------------------------------------------!
|
||||
|
||||
subroutine ions_nose_init( tempw_ , fnosep_ , nhpcl_ , nhptyp_ , ndega_ )
|
||||
use constants, only: factem, pi, terahertz
|
||||
use control_flags, only: tnosep
|
||||
use ions_base, only: ndofp, tions_base_init, nsp, nat, na
|
||||
real(KIND=dbl), intent(in) :: tempw_ , fnosep_(:)
|
||||
integer, intent(in) :: nhpcl_ , nhptyp_ , ndega_
|
||||
integer :: nsvar, i, j, iat, is, ia
|
||||
|
||||
IF( .NOT. tions_base_init ) &
|
||||
CALL errore(' ions_nose_init ', ' you should call ions_base_init first ', 1 )
|
||||
tempw = tempw_
|
||||
fnosep = 0.0d0
|
||||
nhpcl = MAX( nhpcl_ , 1 )
|
||||
nhpdim = 1
|
||||
nhpend = 0
|
||||
atm2nhp(1:nat) = 1
|
||||
nhptyp = 0
|
||||
if (abs(nhptyp_).eq.1) then
|
||||
nhptyp = 1
|
||||
if (nhptyp_.gt.0) nhpend = 1
|
||||
nhpdim = nsp
|
||||
iat = 0
|
||||
do is=1,nsp
|
||||
do ia=1,na(is)
|
||||
iat = iat+1
|
||||
atm2nhp(iat) = is
|
||||
enddo
|
||||
enddo
|
||||
elseif (abs(nhptyp_).eq.2) then
|
||||
nhptyp = 2
|
||||
if (nhptyp_.gt.0) nhpend = 1
|
||||
nhpdim = nat
|
||||
do i=1,nat
|
||||
atm2nhp(i) = i
|
||||
enddo
|
||||
endif
|
||||
! Add one more chain on top if needed
|
||||
nhpdim = nhpdim + nhpend
|
||||
|
||||
IF( nhpcl > nhclm ) &
|
||||
CALL errore(' ions_nose_init ', ' nhpcl out of range ', nhpcl )
|
||||
call ions_nose_allocate()
|
||||
! Setup Nose-Hoover chain masses
|
||||
!
|
||||
if ( ndega_ > 0 ) then
|
||||
ndega = ndega_
|
||||
else if ( ndega_ < 0 ) then
|
||||
ndega = ndofp - ( - ndega_ )
|
||||
else
|
||||
ndega = ndofp
|
||||
endif
|
||||
|
||||
! there is no need to initialize any Nose variables except for nhpcl
|
||||
! and ndega if the thermostat is not used
|
||||
!
|
||||
|
||||
IF( tnosep ) THEN
|
||||
|
||||
IF( nhpcl > SIZE( fnosep_ ) ) &
|
||||
CALL errore(' ions_nose_init ', ' fnosep size too small ', nhpcl )
|
||||
|
||||
! count the number of atoms per thermostat and set the value
|
||||
anum2nhp = 0
|
||||
iat = 0
|
||||
do is=1,nsp
|
||||
do ia=1,na(is)
|
||||
iat = iat+1
|
||||
anum2nhp(atm2nhp(iat)) = anum2nhp(atm2nhp(iat)) + 3
|
||||
enddo
|
||||
enddo
|
||||
if (nhpend.eq.1) anum2nhp(nhpdim) = nhpdim - 1
|
||||
! set gkbt2nhp for each thermostat
|
||||
do is=1,nhpdim
|
||||
gkbt2nhp(is) = DBLE(anum2nhp(is)) * tempw / factem
|
||||
enddo
|
||||
! scale the target energy by some factor convering 3*nat to ndega
|
||||
if (nhpdim.gt.1) then
|
||||
do is=1,(nhpdim-nhpend)
|
||||
gkbt2nhp(is) = gkbt2nhp(is)*DBLE(ndega)/DBLE(3*nat)
|
||||
enddo
|
||||
endif
|
||||
!
|
||||
gkbt = DBLE( ndega ) * tempw / factem
|
||||
if (nhpdim.eq.1) gkbt2nhp(1) = gkbt
|
||||
kbt = tempw / factem
|
||||
|
||||
fnosep(1) = fnosep_ (1)
|
||||
if( fnosep(1) > 0.0d0 ) then
|
||||
qnp_(1) = 2.d0 * gkbt / ( fnosep(1) * ( 2.d0 * pi ) * terahertz )**2
|
||||
end if
|
||||
|
||||
if ( nhpcl > 1 ) then
|
||||
do i = 2, nhpcl
|
||||
fnosep(i) = fnosep_ (i)
|
||||
if( fnosep(i) > 0.0d0 ) then
|
||||
qnp_(i) = 2.d0 * tempw / factem / ( fnosep(i) * ( 2.d0 * pi ) * terahertz )**2
|
||||
else
|
||||
qnp_(i) = qnp_(1) / dble(ndega)
|
||||
endif
|
||||
enddo
|
||||
endif
|
||||
! set the NH masses for all the chains
|
||||
do j=1,nhpdim
|
||||
qnp((j-1)*nhpcl+1) = qnp_(1)*gkbt2nhp(j)/gkbt
|
||||
If (nhpcl > 1) then
|
||||
do i=2,nhpcl
|
||||
qnp((j-1)*nhpcl+i) = qnp_(i)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
END IF
|
||||
|
||||
|
||||
! WRITE( stdout,100)
|
||||
! WRITE( stdout,110) QNOSEP,TEMPW
|
||||
! WRITE( stdout,120) GLIB
|
||||
! WRITE( stdout,130) NSVAR
|
||||
100 FORMAT(//' * Temperature control of ions with nose thermostat'/)
|
||||
110 FORMAT(3X,'nose mass:',F12.4,' temperature (K):',F12.4)
|
||||
120 FORMAT(3X,'ionic degrees of freedom: ',F5.0)
|
||||
130 FORMAT(3X,'time steps per nose oscillation: ',I5,//)
|
||||
|
||||
return
|
||||
end subroutine ions_nose_init
|
||||
|
||||
|
||||
SUBROUTINE ions_nose_allocate()
|
||||
IMPLICIT NONE
|
||||
allocate(vnhp(nhpcl*nhpdim))
|
||||
vnhp = 0.0d0
|
||||
allocate(xnhp0(nhpcl*nhpdim))
|
||||
xnhp0 = 0.0d0
|
||||
allocate(xnhpm(nhpcl*nhpdim))
|
||||
xnhpm = 0.0d0
|
||||
allocate(xnhpp(nhpcl*nhpdim))
|
||||
xnhpp = 0.0d0
|
||||
allocate(ekin2nhp(nhpdim))
|
||||
allocate(gkbt2nhp(nhpdim))
|
||||
allocate(anum2nhp(nhpdim))
|
||||
allocate(qnp(nhpcl*nhpdim))
|
||||
qnp = 0.0d0
|
||||
allocate(qnp_(nhpcl))
|
||||
qnp_ = 0.0d0
|
||||
RETURN
|
||||
END SUBROUTINE ions_nose_allocate
|
||||
|
||||
SUBROUTINE ions_nose_deallocate()
|
||||
IMPLICIT NONE
|
||||
IF( ALLOCATED( vnhp ) ) deallocate(vnhp )
|
||||
IF( ALLOCATED( xnhp0 ) ) deallocate(xnhp0 )
|
||||
IF( ALLOCATED( xnhpm ) ) deallocate(xnhpm )
|
||||
IF( ALLOCATED( xnhpp ) ) deallocate(xnhpp )
|
||||
IF( ALLOCATED( ekin2nhp ) ) deallocate(ekin2nhp )
|
||||
IF( ALLOCATED( gkbt2nhp ) ) deallocate(gkbt2nhp )
|
||||
IF( ALLOCATED( anum2nhp ) ) deallocate(anum2nhp )
|
||||
IF( ALLOCATED( qnp ) ) deallocate(qnp )
|
||||
IF( ALLOCATED( qnp_ ) ) deallocate(qnp_ )
|
||||
RETURN
|
||||
END SUBROUTINE ions_nose_deallocate
|
||||
|
||||
|
||||
SUBROUTINE ions_nose_info()
|
||||
|
||||
use constants, only: factem, terahertz, pi
|
||||
use time_step, only: delt
|
||||
USE io_global, ONLY: stdout
|
||||
USE control_flags, ONLY: tnosep
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
INTEGER :: nsvar, i, j
|
||||
REAL(dbl) :: wnosep
|
||||
|
||||
IF( tnosep ) THEN
|
||||
!
|
||||
IF( fnosep(1) <= 0.D0) &
|
||||
CALL errore(' ions_nose_info ', ' fnosep less than zero ', 1)
|
||||
IF( delt <= 0.D0) &
|
||||
CALL errore(' ions_nose_info ', ' delt less than zero ', 1)
|
||||
|
||||
wnosep = fnosep(1) * ( 2.d0 * pi ) * terahertz
|
||||
nsvar = ( 2.d0 * pi ) / ( wnosep * delt )
|
||||
|
||||
WRITE( stdout,563) tempw, nhpcl, ndega, nsvar
|
||||
WRITE( stdout,564) (fnosep(i),i=1,nhpcl)
|
||||
WRITE( stdout,565) nhptyp, (nhpdim-nhpend), nhpend , &
|
||||
(anum2nhp(j),j=1,nhpdim)
|
||||
do j=1,nhpdim
|
||||
WRITE( stdout,566) j,(qnp((j-1)*nhpcl+i),i=1,nhpcl)
|
||||
enddo
|
||||
END IF
|
||||
|
||||
563 format( //, &
|
||||
& 3X,'ion dynamics with nose` temperature control:', /, &
|
||||
& 3X,'temperature required = ', f10.5, ' (kelvin) ', /, &
|
||||
& 3X,'NH chain length = ', i3, /, &
|
||||
& 3X,'active degrees of freedom = ', i3, /, &
|
||||
& 3X,'time steps per nose osc. = ', i5 )
|
||||
564 format( //, &
|
||||
& 3X,'nose` frequency(es) = ', 20(1X,f10.3) )
|
||||
! 565 format( //, &
|
||||
! & 3X,'nose` mass(es) = ', 20(1X,f10.3), // )
|
||||
565 FORMAT( //, &
|
||||
& 3X,'the requested type of NH chains is ',I5, /, &
|
||||
& 3X,'total number of thermostats used ',I5,1X,I1, /, &
|
||||
& 3X,'ionic degrees of freedom for each chain ',20(1X,I3))
|
||||
566 format( //, &
|
||||
& 3X,'nose` mass(es) for chain ',i4,' = ', 20(1X,f10.3))
|
||||
RETURN
|
||||
END SUBROUTINE ions_nose_info
|
||||
|
||||
|
||||
|
||||
subroutine ions_nosevel( vnhp, xnhp0, xnhpm, delt, nhpcl, nhpdim )
|
||||
implicit none
|
||||
real(KIND=dbl), intent(inout) :: vnhp(nhpcl,nhpdim)
|
||||
real(KIND=dbl), intent(in) :: xnhp0(nhpcl,nhpdim), xnhpm(nhpcl,nhpdim), delt
|
||||
integer, intent(in) :: nhpcl, nhpdim
|
||||
integer :: i,j
|
||||
do j=1,nhpdim
|
||||
do i=1,nhpcl
|
||||
vnhp(i,j)=2.*(xnhp0(i,j)-xnhpm(i,j))/delt-vnhp(i,j)
|
||||
enddo
|
||||
enddo
|
||||
!
|
||||
! this is equivalent to:
|
||||
! velocity = ( 3.D0 * xnos0(1) - 4.D0 * xnosm(1) + xnos2m(1) ) / ( 2.0d0 * delt )
|
||||
! but we do not need variables at time t-2dt ( xnos2m )
|
||||
!
|
||||
return
|
||||
end subroutine ions_nosevel
|
||||
|
||||
|
||||
subroutine ions_noseupd( xnhpp, xnhp0, xnhpm, delt, qnp, ekin2nhp, gkbt2nhp, vnhp, kbt, nhpcl, nhpdim, nhpend )
|
||||
implicit none
|
||||
real(KIND=dbl), intent(out) :: xnhpp(nhpcl,nhpdim), vnhp(nhpcl,nhpdim)
|
||||
real(KIND=dbl), intent(in) :: xnhp0(nhpcl,nhpdim), xnhpm(nhpcl,nhpdim), delt, qnp(nhpcl,nhpdim), gkbt2nhp(:), kbt
|
||||
real(KIND=dbl), intent(inout) :: ekin2nhp(:)
|
||||
integer, intent(in) :: nhpcl, nhpdim, nhpend
|
||||
integer :: i, j
|
||||
real(KIND=dbl) :: dt2, zetfrc, vp1dlt, ekinend, vp1dend
|
||||
|
||||
|
||||
ekinend = 0.0d0
|
||||
vp1dend = 0.0d0
|
||||
if (nhpend.eq.1) vp1dend = 0.5d0*delt*vnhp(1,nhpdim)
|
||||
dt2 = delt**2
|
||||
do j=1,nhpdim
|
||||
zetfrc = dt2*(2.0d0*ekin2nhp(j)-gkbt2nhp(j))
|
||||
If (nhpcl.gt.1) then
|
||||
do i=1,(nhpcl-1)
|
||||
vp1dlt = 0.5d0*delt*vnhp(i+1,j)
|
||||
xnhpp(i,j)=(2.d0*xnhp0(i,j)-(1.d0-vp1dlt)*xnhpm(i,j)+zetfrc/qnp(i,j))&
|
||||
& /(1.d0+vp1dlt)
|
||||
! xnhpp(i,j)=(4.d0*xnhp0(i,j)-(2.d0-delt*vnhp(i+1,j))*xnhpm(i,j)+2.0d0*dt2*zetfrc/qnp(i,j))&
|
||||
! & /(2.d0+delt*vnhp(i+1,j))
|
||||
vnhp(i,j) =(xnhpp(i,j)-xnhpm(i,j))/( 2.0d0 * delt )
|
||||
zetfrc = dt2*(qnp(i,j)*vnhp(i,j)**2-kbt)
|
||||
enddo
|
||||
endif
|
||||
! Last variable
|
||||
i = nhpcl
|
||||
if (nhpend.eq.0) then
|
||||
xnhpp(i,j)=2.d0*xnhp0(i,j)-xnhpm(i,j)+ zetfrc / qnp(i,j)
|
||||
vnhp(i,j) =(xnhpp(i,j)-xnhpm(i,j))/( 2.0d0 * delt )
|
||||
elseif (nhpend.eq.1) then
|
||||
xnhpp(i,j)=(2.d0*xnhp0(i,j)-(1.d0-vp1dend)*xnhpm(i,j)+zetfrc/qnp(i,j))&
|
||||
& /(1.d0+vp1dend)
|
||||
vnhp(i,j) =(xnhpp(i,j)-xnhpm(i,j))/( 2.0d0 * delt )
|
||||
ekinend = ekinend + (qnp(i,j)*vnhp(i,j)**2)
|
||||
if (j.eq.(nhpdim-nhpend)) then
|
||||
ekin2nhp(nhpdim) = 0.5d0*ekinend
|
||||
vp1dend = 0.0d0
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
! Update velocities
|
||||
! do i=1,nhpcl
|
||||
! vnhp(i) =(xnhpp(i)-xnhpm(i))/( 2.0d0 * delt )
|
||||
! end do
|
||||
! These are the original expressions from cpr.f90
|
||||
! xnhpp(1)=2.*xnhp0(1)-xnhpm(1)+2.*( delt**2 / qnp(1) )*(ekinpr-gkbt/2.)
|
||||
! vnhp(1) =(xnhpp(1)-xnhpm(1))/( 2.0d0 * delt )
|
||||
return
|
||||
end subroutine ions_noseupd
|
||||
|
||||
|
||||
|
||||
real(KIND=dbl) function ions_nose_nrg( xnhp0, vnhp, qnp, gkbt2nhp, kbt, nhpcl, nhpdim )
|
||||
implicit none
|
||||
integer :: nhpcl, nhpdim
|
||||
real(KIND=dbl) :: gkbt2nhp(:), qnp(nhpcl,nhpdim),vnhp(nhpcl,nhpdim),xnhp0(nhpcl,nhpdim),kbt
|
||||
integer :: i,j
|
||||
real(KIND=dbl) :: stmp
|
||||
!
|
||||
stmp = 0.0d0
|
||||
do j=1,nhpdim
|
||||
stmp = stmp + 0.5d0 * qnp(1,j) * vnhp(1,j) * vnhp(1,j) + gkbt2nhp(j) * xnhp0(1,j)
|
||||
if (nhpcl > 1) then
|
||||
do i=2,nhpcl
|
||||
stmp=stmp+0.5*qnp(i,j)*vnhp(i,j)*vnhp(i,j) + kbt*xnhp0(i,j)
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
ions_nose_nrg = stmp
|
||||
return
|
||||
end function ions_nose_nrg
|
||||
|
||||
|
||||
|
||||
subroutine ions_nose_shiftvar( xnhpp, xnhp0, xnhpm )
|
||||
! shift values of nose variables to start a new step
|
||||
implicit none
|
||||
real(KIND=dbl), intent(inout) :: xnhpm(:), xnhp0(:)
|
||||
real(KIND=dbl), intent(in) :: xnhpp(:)
|
||||
!
|
||||
xnhpm = xnhp0
|
||||
xnhp0 = xnhpp
|
||||
!
|
||||
return
|
||||
end subroutine ions_nose_shiftvar
|
||||
|
||||
!------------------------------------------------------------------------------!
|
||||
END MODULE ions_nose
|
||||
!------------------------------------------------------------------------------!
|
Loading…
Reference in New Issue