scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

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


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1998 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
kkudin 2005-07-05 21:02:48 +00:00
parent b84ae73507
commit 930619e2b2
11 changed files with 231 additions and 103 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
CPV/cp_restart.f90
View File

@ -54,6 +54,7 @@ MODULE cp_restart
USE ions_base, ONLY: nsp, nat, na, atm, zv, pmass, amass, iforce
USE funct, ONLY: dft
USE mp, ONLY: mp_sum
USE parameters, ONLY: nhclm
IMPLICIT NONE
INTEGER, INTENT(IN) :: ndw !
@ -271,8 +272,8 @@ MODULE cp_restart
!
call iotk_write_begin(iunpun,"IONS_NOSE")
call iotk_write_dat (iunpun, "nhpcl", nhpcl)
call iotk_write_dat (iunpun, "xnhp", xnhp0 )
call iotk_write_dat (iunpun, "vnhp", vnhp)
call iotk_write_dat (iunpun, "xnhp", xnhp0(1:nhclm) )
call iotk_write_dat (iunpun, "vnhp", vnhp(1:nhclm) )
call iotk_write_end(iunpun,"IONS_NOSE")
!
call iotk_write_dat (iunpun, "ekincm", ekincm)
@ -305,7 +306,7 @@ MODULE cp_restart
!
call iotk_write_begin(iunpun,"IONS_NOSE")
call iotk_write_dat (iunpun, "nhpcl", nhpcl)
call iotk_write_dat (iunpun, "xnhp", xnhpm )
call iotk_write_dat (iunpun, "xnhp", xnhpm(1:nhclm) )
! call iotk_write_dat (iunpun, "vnhp", vnhp)
call iotk_write_end(iunpun,"IONS_NOSE")
!
@ -569,6 +570,7 @@ MODULE cp_restart
USE ions_base, ONLY: nsp, nat, na, atm, zv, pmass
USE reciprocal_vectors, ONLY: ngwt, ngw, ig_l2g, mill_l
USE mp, ONLY: mp_sum, mp_bcast
USE parameters, ONLY: nhclm
IMPLICIT NONE
INTEGER, INTENT(IN) :: ndr ! I/O unit number
@ -689,8 +691,8 @@ MODULE cp_restart
!
call iotk_scan_begin(iunpun,"IONS_NOSE")
call iotk_scan_dat (iunpun, "nhpcl", nhpcl_ )
call iotk_scan_dat (iunpun, "xnhp", xnhp0 )
call iotk_scan_dat (iunpun, "vnhp", vnhp)
call iotk_scan_dat (iunpun, "xnhp", xnhp0(1:nhclm) )
call iotk_scan_dat (iunpun, "vnhp", vnhp(1:nhclm) )
call iotk_scan_end(iunpun,"IONS_NOSE")
!
call iotk_scan_dat (iunpun, "ekincm", ekincm)
@ -730,7 +732,7 @@ MODULE cp_restart
!
call iotk_scan_begin(iunpun,"IONS_NOSE")
call iotk_scan_dat (iunpun, "nhpcl", nhpcl_ )
call iotk_scan_dat (iunpun, "xnhp", xnhpm )
call iotk_scan_dat (iunpun, "xnhp", xnhpm(1:nhclm) )
call iotk_scan_end(iunpun,"IONS_NOSE")
!
call iotk_scan_begin(iunpun,"ELECTRONS_NOSE")

42
CPV/cpr.f90
View File

@ -126,10 +126,13 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
USE cp_electronic_mass, ONLY : emass, emass_cutoff, emass_precond
USE ions_positions, ONLY : tau0, taum, taup, taus, tausm, tausp, &
vels, velsm, velsp, ions_hmove, ions_move
USE ions_nose, ONLY : gkbt, kbt, ndega, nhpcl, qnp, vnhp, &
xnhp0, xnhpm, xnhpp, ions_nosevel, &
ions_noseupd, tempw, ions_nose_nrg, &
ions_nose_shiftvar
USE ions_nose, ONLY : gkbt, kbt, ndega, nhpcl, nhpdim, qnp, &
vnhp, xnhp0, xnhpm, xnhpp, atm2nhp, &
ions_nosevel, ions_noseupd, &
ions_nose_allocate, ions_nose_deallocate,&
tempw, ions_nose_nrg, &
ions_nose_shiftvar, gkbt2nhp, ekin2nhp, &
anum2nhp
USE electrons_nose, ONLY : qne, ekincw, xnhe0, xnhep, xnhem, &
vnhe, electrons_nose_nrg, &
electrons_nose_shiftvar, &
@ -435,7 +438,7 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
!
IF ( .NOT. tsde ) fccc = 1.D0 / ( 1.D0 + frice )
!
IF ( tnosep ) CALL ions_nosevel( vnhp, xnhp0, xnhpm, delt, nhpcl )
IF ( tnosep ) CALL ions_nosevel( vnhp, xnhp0, xnhpm, delt, nhpcl, nhpdim )
!
IF ( tnosee ) THEN
!
@ -525,7 +528,8 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
!
CALL ions_move( tausp, taus, tausm, iforce, pmass, fion, &
ainv, delt, na, nsp, fricp, hgamma, vels, &
tsdp, tnosep, fionm, vnhp, velsp, velsm )
tsdp, tnosep, fionm, vnhp, velsp, velsm, &
nhpcl, nhpdim, atm2nhp )
!
IF ( lconstrain ) THEN
!
@ -641,7 +645,8 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
! ... ionic temperature
!
IF ( tfor ) CALL ions_temp( tempp, temps, ekinpr, vels, &
na, nsp, hold, pmass, ndega )
na, nsp, hold, pmass, ndega, &
nhpdim, atm2nhp, ekin2nhp )
!
! ... fake electronic kinetic energy
!
@ -672,7 +677,7 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
! ... udating nose-hoover friction variables
!
IF ( tnosep ) CALL ions_noseupd( xnhpp, xnhp0, xnhpm, delt, &
qnp, ekinpr, gkbt, vnhp, kbt, nhpcl )
qnp, ekin2nhp, gkbt2nhp, vnhp, kbt, nhpcl, nhpdim )
!
IF ( tnosee ) CALL electrons_noseupd( xnhep, xnhe0, xnhem, &
delt, qne, ekinc, ekincw, vnhe )
@ -751,7 +756,7 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
! ... add energies of thermostats
!
IF ( tnosep ) econt = econt + &
ions_nose_nrg( xnhp0, vnhp, qnp, gkbt, kbt, nhpcl )
ions_nose_nrg( xnhp0, vnhp, qnp, gkbt2nhp, kbt, nhpcl, nhpdim )
IF ( tnosee ) econt = econt + &
electrons_nose_nrg( xnhe0, vnhe, qne, ekincw )
IF ( tnoseh ) econt = econt + &
@ -806,12 +811,12 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
!
WRITE( stdout,11)
!
CALL printout_pos( stdout, nfi, tau0, nat, tps )
! CALL printout_pos( stdout, nfi, tau0, nat, tps )
CALL printout_pos( 35 , nfi, tau0, nat, tps )
!
! ... write out a standard XYZ file in angstroms
!
! CALL print_pos_in( stdout, nfi, tau0 , nat, tps, ityp, atm,ind_bck,one)
CALL print_pos_in( stdout, nfi, tau0 , nat, tps, ityp, atm,ind_bck,one)
! WRITE( 35, '(I5)') nat
! CALL print_pos_in( 35 , nfi, tau0 , nat, tps, ityp, atm,ind_bck,toang)
@ -834,16 +839,16 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
!
WRITE( stdout, 12 )
!
CALL printout_pos( stdout, nfi, tauw, nat, tps )
! CALL printout_pos( stdout, nfi, tauw, nat, tps )
CALL printout_pos( 34 , nfi, tauw, nat, tps )
! CALL print_pos_in( stdout, nfi, tauw , nat, tps, ityp, atm,ind_bck,one)
CALL print_pos_in( stdout, nfi, tauw , nat, tps, ityp, atm,ind_bck,one)
! CALL print_pos_in( 34 , nfi, tauw , nat, tps, ityp, atm,ind_bck,one)
!
WRITE( stdout, 13 )
!
CALL printout_pos( stdout, nfi, fion, nat, tps )
! CALL printout_pos( stdout, nfi, fion, nat, tps )
CALL printout_pos( 37 , nfi, fion, nat, tps )
! CALL print_pos_in( stdout, nfi, fion , nat, tps, ityp, atm,ind_bck,one)
CALL print_pos_in( stdout, nfi, fion , nat, tps, ityp, atm,ind_bck,one)
! CALL print_pos_in( 37 , nfi, fion , nat, tps, ityp, atm,ind_bck,one)
!
DEALLOCATE( tauw )
@ -1101,6 +1106,13 @@ SUBROUTINE cprmain( tau, fion_out, etot_out )
!
1977 FORMAT(5X,//'====================== end cprvan ======================',//)
!
! by Kostya
! Something is fishy here, when deallocate_modules_var is called
! IFC 8.0 and 8.1 spit out
!*** glibc detected *** free(): invalid next size (fast): 0x18b46af8 ***
!forrtl: error (76): IOT trap signal
! I could not find what is wrong ...
call ions_nose_deallocate()
CALL deallocate_modules_var()
!
RETURN

6
CPV/input.f90
View File

@ -752,8 +752,8 @@ MODULE input
noptical, boptical, k_points, nkstot, nk1, nk2, nk3, k1, k2, k3, &
xk, wk, occupations, n_inner, fermi_energy, rotmass, occmass, &
rotation_damping, occupation_damping, occupation_dynamics, &
rotation_dynamics, degauss, smearing, nhpcl, ndega, cell_units, &
restart_mode
rotation_dynamics, degauss, smearing, nhpcl, nhptyp, ndega, &
cell_units, restart_mode
USE input_parameters, ONLY: diis_achmix, diis_ethr, diis_wthr, diis_delt, &
diis_nreset, diis_temp, diis_nrot, diis_maxstep, diis_fthr, &
@ -880,7 +880,7 @@ MODULE input
CALL cell_nose_init( temph, fnoseh )
CALL ions_nose_init( tempw, fnosep, nhpcl, ndega, nat )
CALL ions_nose_init( tempw, fnosep, nhpcl, nhptyp, ndega )
CALL electrons_nose_init( ekincw , fnosee )

220
CPV/ions_nose.f90
View File

@ -10,16 +10,35 @@
!------------------------------------------------------------------------------!
USE kinds, ONLY: dbl
USE parameters, ONLY: nhclm
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
! 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
INTEGER :: atm2nhp(natx)
INTEGER, ALLOCATABLE :: anum2nhp(:)
REAL(dbl), ALLOCATABLE :: vnhp(:), xnhp0(:), xnhpm(:), xnhpp(:), &
ekin2nhp(:), gkbt2nhp(:), qnp(:), qnp_(:)
INTEGER :: nhpcl, ndega
REAL(dbl) :: vnhp( nhclm ) = 0.0d0
REAL(dbl) :: xnhp0( nhclm ) = 0.0d0
REAL(dbl) :: xnhpm( nhclm ) = 0.0d0
REAL(dbl) :: xnhpp( nhclm ) = 0.0d0
REAL(dbl) :: qnp( nhclm ) = 0.0d0
REAL(dbl) :: gkbt = 0.0d0
REAL(dbl) :: kbt = 0.0d0
REAL(dbl) :: tempw = 0.0d0
@ -29,33 +48,42 @@
CONTAINS
!------------------------------------------------------------------------------!
subroutine ions_nose_init( tempw_ , fnosep_ , nhpcl_ , ndega_ , nat )
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
implicit none
use ions_base, only: ndofp, tions_base_init, nsp, nat, na
real(KIND=dbl), intent(in) :: tempw_ , fnosep_(:)
integer, intent(in) :: nhpcl_ , ndega_
integer, intent(in) :: nat
integer :: nsvar, i
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 )
vnhp = 0.0d0
xnhp0 = 0.0d0
xnhpm = 0.0d0
xnhpp = 0.0d0
tempw = tempw_
fnosep = 0.0d0
qnp = 0.0d0
nhpcl = MAX( nhpcl_ , 1 )
nhpdim = 1
atm2nhp(1:nat) = 1
nhptyp = 0
if (nhptyp_.eq.1) then
nhptyp = 1
nhpdim = nsp
iat = 0
do is=1,nsp
do ia=1,na(is)
iat = iat+1
atm2nhp(iat) = is
enddo
enddo
elseif (nhptyp_.eq.2) then
nhptyp = 2
nhpdim = nat
do i=1,nat
atm2nhp(i) = i
enddo
endif
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
@ -75,26 +103,49 @@
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
! set gkbt2nhp for each thermostat
do is=1,nhpdim
gkbt2nhp(is) = DBLE(anum2nhp(is)) * tempw / factem
enddo
!
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
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
qnp_(i) = 2.d0 * tempw / factem / ( fnosep(i) * ( 2.d0 * pi ) * terahertz )**2
else
qnp(i) = qnp(1) / dble(ndega)
qnp_(i) = qnp_(1) / dble(ndega)
endif
enddo
endif
END IF
! 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)
@ -110,6 +161,40 @@
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( qnp ) ) deallocate(qnp )
IF( ALLOCATED( qnp_ ) ) deallocate(qnp_ )
RETURN
END SUBROUTINE ions_nose_deallocate
SUBROUTINE ions_nose_info()
use constants, only: factem, terahertz, pi
@ -119,7 +204,7 @@
IMPLICIT NONE
INTEGER :: nsvar, i
INTEGER :: nsvar, i, j
REAL(dbl) :: wnosep
IF( tnosep ) THEN
@ -134,8 +219,11 @@
WRITE( stdout,563) tempw, nhpcl, ndega, nsvar
WRITE( stdout,564) (fnosep(i),i=1,nhpcl)
WRITE( stdout,565) (qnp(i),i=1,nhpcl)
END IF
WRITE( stdout,565) nhptyp, nhpdim, (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:', /, &
@ -145,23 +233,29 @@
& 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,'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, /, &
& 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 )
subroutine ions_nosevel( vnhp, xnhp0, xnhpm, delt, nhpcl, nhpdim )
implicit none
real(KIND=dbl), intent(inout) :: vnhp(:)
real(KIND=dbl), intent(in) :: xnhp0(:), xnhpm(:), delt
integer, intent(in) :: nhpcl
integer :: i
do i=1,nhpcl
vnhp(i)=2.*(xnhp0(i)-xnhpm(i))/delt-vnhp(i)
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:
@ -172,28 +266,31 @@
end subroutine ions_nosevel
subroutine ions_noseupd( xnhpp, xnhp0, xnhpm, delt, qnp, ekinpr, gkbt, vnhp, kbt, nhpcl )
subroutine ions_noseupd( xnhpp, xnhp0, xnhpm, delt, qnp, ekin2nhp, gkbt2nhp, vnhp, kbt, nhpcl, nhpdim )
implicit none
real(KIND=dbl), intent(out) :: xnhpp(:), vnhp(:)
real(KIND=dbl), intent(in) :: xnhp0(:), xnhpm(:), delt, qnp(:), ekinpr, gkbt, kbt
integer, intent(in) :: nhpcl
integer :: i
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), ekin2nhp(:), gkbt2nhp(:), kbt
integer, intent(in) :: nhpcl, nhpdim
integer :: i, j
real(KIND=dbl) :: dt2, zetfrc
zetfrc = 2.0d0*ekinpr-gkbt
dt2 = delt**2
do j=1,nhpdim
zetfrc = 2.0d0*ekin2nhp(j)-gkbt2nhp(j)
If (nhpcl.gt.1) then
do i=1,(nhpcl-1)
xnhpp(i)=(4.d0*xnhp0(i)-(2.d0-delt*vnhp(i+1))*xnhpm(i)+2.0d0*dt2*zetfrc/qnp(i))&
& /(2.d0+delt*vnhp(i+1))
vnhp(i) =(xnhpp(i)-xnhpm(i))/( 2.0d0 * delt )
zetfrc = (qnp(i)*vnhp(i)**2-kbt)
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 = (qnp(i,j)*vnhp(i,j)**2-kbt)
enddo
endif
! Last variable
xnhpp(nhpcl)=2.d0*xnhp0(nhpcl)-xnhpm(nhpcl)+( delt**2 / qnp(nhpcl) )*zetfrc
i = nhpcl
vnhp(i) =(xnhpp(i)-xnhpm(i))/( 2.0d0 * delt )
xnhpp(i,j)=2.d0*xnhp0(i,j)-xnhpm(i,j)+( delt**2 / qnp(i,j) )*zetfrc
vnhp(i,j) =(xnhpp(i,j)-xnhpm(i,j))/( 2.0d0 * delt )
enddo
! Update velocities
! do i=1,nhpcl
! vnhp(i) =(xnhpp(i)-xnhpm(i))/( 2.0d0 * delt )
@ -206,19 +303,22 @@
real(KIND=dbl) function ions_nose_nrg( xnhp0, vnhp, qnp, gkbt, kbt, nhpcl )
real(KIND=dbl) function ions_nose_nrg( xnhp0, vnhp, qnp, gkbt2nhp, kbt, nhpcl, nhpdim )
implicit none
integer :: nhpcl
real(KIND=dbl) :: gkbt,qnp(:),vnhp(:),xnhp0(:),kbt
integer :: i
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.5d0 * qnp(1) * vnhp(1) * vnhp(1) + gkbt * xnhp0(1)
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)*vnhp(i)*vnhp(i) + kbt*xnhp0(i)
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

8
CPV/ions_positions.f90
View File

@ -68,7 +68,7 @@ MODULE ions_positions
!--------------------------------------------------------------------------
SUBROUTINE ions_move( tausp, taus, tausm, iforce, pmass, fion, ainv, &
delt, na, nsp, fricp, hgamma, vels, tsdp, tnosep, &
fionm, vnhp, velsp, velsm )
fionm, vnhp, velsp, velsm, nhpcl, nhpdim, atm2nhp )
!--------------------------------------------------------------------------
!
IMPLICIT NONE
@ -77,11 +77,11 @@ MODULE ions_positions
INTEGER, INTENT(IN) :: iforce(:,:)
REAL(KIND=dbl), INTENT(IN) :: ainv(3,3), delt
REAL(KIND=dbl), INTENT(OUT) :: tausp(:,:)
INTEGER, INTENT(IN) :: na(:), nsp
INTEGER, INTENT(IN) :: na(:), nsp, nhpcl, nhpdim, atm2nhp(:)
REAL(KIND=dbl), INTENT(IN) :: fricp, hgamma(3,3), vels(:,:)
LOGICAL, INTENT(IN) :: tsdp, tnosep
REAL(KIND=dbl), INTENT(INOUT) :: fionm(:,:)
REAL(KIND=dbl), INTENT(IN) :: vnhp(:)
REAL(KIND=dbl), INTENT(IN) :: vnhp(nhpcl,nhpdim)
REAL(KIND=dbl), INTENT(OUT) :: velsp(:,:)
REAL(KIND=dbl), INTENT(IN) :: velsm(:,:)
INTEGER :: is, ia, i, isa
@ -135,7 +135,7 @@ MODULE ions_positions
fionm(i,isa) = ainv(i,1) * fion(1,isa) + &
ainv(i,2) * fion(2,isa) + &
ainv(i,3) * fion(3,isa) - &
vnhp(1) * vels(i,isa) * pmass(is) - &
vnhp(1,atm2nhp(isa)) * vels(i,isa) * pmass(is) - &
pmass(is) * ( hgamma(i,1) * vels(1,isa) + &
hgamma(i,2) * vels(2,isa) + &
hgamma(i,3) * vels(3,isa) )

7
CPV/main.f90
View File

@ -167,7 +167,7 @@
USE smooth_grid_dimensions, ONLY: nr1s, nr2s, nr3s, nr1sx, nr2sx, nr3sx
!
USE ions_nose, ONLY: ions_nose_shiftvar, vnhp, xnhpp, xnhp0, xnhpm, ions_nosevel, &
ions_noseupd, qnp, gkbt, kbt, nhpcl
ions_noseupd, qnp, gkbt, kbt, nhpcl, nhpdim, gkbt2nhp, ekin2nhp
USE wave_init, ONLY: pw_atomic_init
USE electrons_module, ONLY: electron_mass_init, band_init
@ -586,14 +586,15 @@
IF(tfor .AND. doions) THEN
! ... Determines DXNOS/DT dynamically
IF (tnosep) THEN
CALL ions_nosevel( vnhp, xnhp0, xnhpm, delt, 1 )
CALL ions_nosevel( vnhp, xnhp0, xnhpm, delt, 1, 1 )
vnosep = vnhp(1)
END IF
! ... move ionic degrees of freedom
! ... WRITE( stdout,*) '* TSDP *', tsdp
ekinp = moveions(tsdp, thdyn, nfi, atoms_m, atoms_0, atoms_p, ht_m, ht_0, vnosep)
IF (tnosep) THEN
CALL ions_noseupd( xnhpp, xnhp0, xnhpm, delt, qnp, atoms_0%ekint, gkbt, vnhp, kbt, nhpcl )
! below one really should have atoms_0%ekint and NOT ekin2nhp
CALL ions_noseupd( xnhpp, xnhp0, xnhpm, delt, qnp, ekin2nhp, gkbt2nhp, vnhp, kbt, nhpcl, nhpdim )
END IF
END IF

6
CPV/print_out.f90
View File

@ -59,7 +59,7 @@
USE sic_module, ONLY: ind_localisation, pos_localisation, nat_localisation, self_interaction
USE sic_module, ONLY: rad_localisation
USE ions_base, ONLY: ions_temp, cdmi, taui
USE ions_nose, ONLY: ndega, ions_nose_nrg, xnhp0, vnhp, qnp, gkbt, kbt, nhpcl
USE ions_nose, ONLY: ndega, ions_nose_nrg, xnhp0, vnhp, qnp, gkbt, kbt, nhpcl, nhpdim, atm2nhp, ekin2nhp, gkbt2nhp
USE cell_nose, ONLY: cell_nose_nrg, qnh, temph, xnhh0, vnhh
USE cell_base, ONLY: iforceh
USE printout_base, ONLY: printout_base_open, printout_base_close, &
@ -105,7 +105,7 @@
! ... Calculate Ions temperature tempp (in Kelvin )
CALL ions_temp( tempp, temps, ekinpr, atoms%vels, atoms%na, atoms%nsp, ht%hmat, atoms%m, ndega )
CALL ions_temp( tempp, temps, ekinpr, atoms%vels, atoms%na, atoms%nsp, ht%hmat, atoms%m, ndega, nhpdim, atm2nhp, ekin2nhp )
! ... Calculate MSD for each specie, starting from taui positions
@ -133,7 +133,7 @@
endif
IF( tnosep ) THEN
enosep = ions_nose_nrg( xnhp0, vnhp, qnp, gkbt, kbt, nhpcl )
enosep = ions_nose_nrg( xnhp0, vnhp, qnp, gkbt2nhp, kbt, nhpcl, nhpdim )
ELSE
enosep = 0
END IF

9
CPV/smcp.f90
View File

@ -105,7 +105,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
use ions_positions, only: tau0, velsp
use ions_positions, only: ions_hmove, ions_move
use ions_nose, only: gkbt, qnp, ions_nosevel, ions_noseupd, tempw, &
ions_nose_nrg, kbt, nhpcl, ndega
ions_nose_nrg, kbt, nhpcl, ndega, nhpdim, atm2nhp
USE cell_base, ONLY: cell_kinene, cell_move, cell_gamma, cell_hmove
USE cell_nose, ONLY: xnhh0, xnhhm, xnhhp, vnhh, temph, qnh, &
cell_nosevel, cell_noseupd, cell_nose_nrg, cell_nosezero
@ -1174,7 +1174,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
CALL ions_move( rep(sm_k)%tausp, rep(sm_k)%taus, rep(sm_k)%tausm, iforce, pmass, &
rep(sm_k)%fion, ainv, delt, na, nsp, fricp, hgamma, rep(sm_k)%vels, tsdp, &
tnosep, rep(sm_k)%fionm, vnhp(:,sm_k), velsp, rep(sm_k)%velsm )
tnosep, rep(sm_k)%fionm, vnhp(:,sm_k), velsp, rep(sm_k)%velsm, 1, 1, atm2nhp )
!
!cc call cofmass(velsp,rep(sm_k)%cdmvel)
! call cofmass(rep(sm_k)%tausp,cdm)
@ -1328,7 +1328,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
!
if(tnosep)then
CALL ions_noseupd( xnhpp( :, sm_k), xnhp0( :, sm_k), xnhpm( :, sm_k), delt, qnp, &
ekinpr(sm_k), gkbt, vnhp( :, sm_k), kbt, nhpcl )
ekinpr(sm_k), gkbt, vnhp( :, sm_k), kbt, nhpcl, 1 )
endif
if(tnosee)then
call elec_noseupd( xnhep(sm_k), xnhe0(sm_k), xnhem(sm_k), delt, qne, ekinc(sm_k), ekincw, vnhe(sm_k) )
@ -1372,7 +1372,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
econt(sm_k)=econs(sm_k)+ekinc(sm_k)
!
if(tnosep)then
econt(sm_k)=econt(sm_k)+ ions_nose_nrg( xnhp0(:,sm_k), vnhp(:,sm_k), qnp, gkbt, kbt, nhpcl )
econt(sm_k)=econt(sm_k)+ ions_nose_nrg( xnhp0(:,sm_k), vnhp(:,sm_k), qnp, gkbt, kbt, nhpcl, 1 )
endif
if(tnosee)then
econt(sm_k)=econt(sm_k)+ electrons_nose_nrg( xnhe0(sm_k), vnhe(sm_k), qne, ekincw )
@ -1742,6 +1742,7 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out )
IF( ALLOCATED( p_tan )) DEALLOCATE( p_tan )
call ions_nose_deallocate()
CALL deallocate_modules_var()
CALL deallocate_smd_rep()

5
Modules/input_parameters.f90
View File

@ -937,6 +937,9 @@ MODULE input_parameters
! non-zero only with "ion_temperature = 'nose' "
! this defines the length of the Nose-Hoover chain
INTEGER :: nhptyp = 0
! this parameter set the nose hoover thermostat to more than one
INTEGER :: ndega = 0
! this is the parameter to control active degrees of freedom
! used for temperature control and the Nose-Hoover chains
@ -1127,7 +1130,7 @@ MODULE input_parameters
NAMELIST / ions / phase_space, ion_dynamics, ion_radius, ion_damping, &
ion_positions, ion_velocities, ion_temperature, &
tempw, fnosep, nhpcl, ndega, tranp, amprp, greasp, &
tempw, fnosep, nhpcl, nhptyp, ndega, tranp, amprp, greasp, &
tolp, ion_nstepe, ion_maxstep, upscale, delta_t, &
pot_extrapolation, wfc_extrapolation, nraise, &
num_of_images, CI_scheme, opt_scheme, use_masses, &

15
Modules/ions_base.f90
View File

@ -578,21 +578,23 @@
!------------------------------------------------------------------------------!
subroutine ions_temp( tempp, temps, ekinpr, vels, na, nsp, h, pmass, ndega )
subroutine ions_temp( tempp, temps, ekinpr, vels, na, nsp, h, pmass, ndega, nhpdim, atm2nhp, ekin2nhp )
use constants, only: factem
implicit none
real( kind=8 ), intent(out) :: ekinpr, tempp
real( kind=8 ), intent(out) :: temps(:)
real( kind=8 ), intent(out) :: ekin2nhp(:)
real( kind=8 ), intent(in) :: vels(:,:)
real( kind=8 ), intent(in) :: pmass(:)
real( kind=8 ), intent(in) :: h(:,:)
integer, intent(in) :: na(:), nsp, ndega
integer, intent(in) :: na(:), nsp, ndega, nhpdim, atm2nhp(:)
integer :: nat, i, j, is, ia, ii, isa
real( kind=8 ) :: cdmvel(3), eks
real( kind=8 ) :: cdmvel(3), eks, eks1
call ions_cofmass( vels, pmass, na, nsp, cdmvel )
nat = SUM( na(1:nsp) )
ekinpr = 0.0d0
temps( 1:nsp ) = 0.0d0
ekin2nhp(1:nhpdim) = 0.0d0
do i=1,3
do j=1,3
do ii=1,3
@ -601,7 +603,9 @@
eks = 0.0d0
do ia=1,na(is)
isa = isa + 1
eks=eks+pmass(is)*h(j,i)*(vels(i,isa)-cdmvel(i))*h(j,ii)*(vels(ii,isa)-cdmvel(ii))
eks1 = pmass(is)*h(j,i)*(vels(i,isa)-cdmvel(i))*h(j,ii)*(vels(ii,isa)-cdmvel(ii))
eks=eks+eks1
ekin2nhp(atm2nhp(isa)) = ekin2nhp(atm2nhp(isa)) + eks1
end do
ekinpr = ekinpr + eks
temps(is) = temps(is) + eks
@ -609,6 +613,9 @@
end do
end do
end do
do is = 1, nhpdim
ekin2nhp(is) = ekin2nhp(is) * 0.5d0
enddo
do is = 1, nsp
temps( is ) = temps( is ) * 0.5d0
temps( is ) = temps( is ) * factem / ( 1.5d0 * na(is) )

2
Modules/read_namelists.f90
View File

@ -338,6 +338,7 @@ MODULE read_namelists_module
fnosep = -1.0D0
fnosep(1) = 1.0D0
nhpcl = 0
nhptyp = 0
ndega = 0
tranp = .FALSE.
amprp = 0.D0
@ -797,6 +798,7 @@ MODULE read_namelists_module
CALL mp_bcast( tempw, ionode_id )
CALL mp_bcast( fnosep, ionode_id )
CALL mp_bcast( nhpcl, ionode_id )
CALL mp_bcast( nhptyp, ionode_id )
CALL mp_bcast( ndega, ionode_id )
CALL mp_bcast( tranp, ionode_id )
CALL mp_bcast( amprp, ionode_id )