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

wfc-extrapolation extended to all "relax" algorithms. 

Molecular Dynamics based algorithm partially rewritten:
both standard and damped MD are performed with the velocity Verlet
scheme (with or without constrains).
Renata's subroutines are used only in the framework of variable
cell.
constrain.f90 file is no longer needed: cnstrains are set in the
input file (see CONSTRAINTS CARD) with the same input format used
in FPMD. An arbitrary number of constrains can be set.
In the case of constrained relaxation the damped MD
algorithm is used instead of BFGS.
When restart_mode = "from_scratch" many reastr files are removed
from the scratch directory.
Several other modifications here and there. (C.S.)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@675 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
sbraccia 2004-02-26 11:50:36 +00:00
parent 5eb515bbe8
commit 0224c2a9ad
13 changed files with 1346 additions and 1292 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Modules/control_flags.f90
View File

@ -213,8 +213,10 @@
lmd, &! if .TRUE. the calculation is a dynamics
lneb, &! if .TRUE. the calculation is neb
lphonon, &! if .TRUE. the calculation is phonon
lconstrain, &! if .TRUE. the calculation is constraint
ldamped, &! if .TRUE. the calculation is a damped dynamics
conv_elec, &! if .TRUE. electron convergence has been reached
conv_ions, &! if .TRUE. ionic convergence has been reached
conv_ions, &! if .TRUE. ionic convergence has been reached
nosym, &! if .TRUE. no symmetry is used
noinv = .FALSE., &! if .TRUE. eliminates inversion symmetry
diis_wfc_keep, &! if .TRUE. keeps old wfc for starting

8
Modules/input_parameters.f90
View File

@ -726,7 +726,8 @@ MODULE input_parameters
!
CHARACTER(LEN=80) :: ion_dynamics = 'none'
! ion_dynamics = 'sd' | 'cg' | 'damp' | 'verlet' | 'bfgs' | 'old-bfgs' | 'none'*
! ion_dynamics = 'sd' | 'cg' | 'damp' | 'verlet' | 'bfgs' |
! 'old-bfgs' | 'none'*
! set how ions shold be moved
! 'none' ions are kept fixed
! 'bfgs' a new BFGS algorithm is used to minimize ionic configuration
@ -737,8 +738,9 @@ MODULE input_parameters
! 'verlet' standard Verlet algorithm is used to propagate ions
CHARACTER(LEN=80) :: ion_dynamics_allowed(10)
DATA ion_dynamics_allowed / 'sd', 'cg', 'damp', 'verlet', 'none', &
'bfgs', 'old-bfgs', 'constrained-bfgs', 'constrained-verlet', 'beeman' /
DATA ion_dynamics_allowed / 'sd', 'cg', 'damp', 'verlet', 'none', &
'bfgs', 'old-bfgs', 'constrained-damp', &
'constrained-verlet', 'beeman' /
REAL(dbl) :: ion_radius(nsx) = 0.5d0
! pseudo-atomic radius of the i-th atomic species

118
Modules/parser.f90
View File

@ -1,55 +1,101 @@
!
! Copyright (C) 2001 Carlo Cavazzoni and PWSCF group
! Copyright (C) 2001-2004 Carlo Cavazzoni and PWSCF group
! 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 parser
!----------------------------------------------------------------------------
MODULE parser
!----------------------------------------------------------------------------
!
USE io_global, ONLY : stdout
USE kinds
!
INTERFACE parser_error
MODULE PROCEDURE p_err_I, p_err_R, p_err_S, p_err_L
END INTERFACE
contains
PURE FUNCTION int_to_char( int )
!
CONTAINS
!
!-----------------------------------------------------------------------
PURE FUNCTION int_to_char( int )
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: int
CHARACTER (LEN=6) :: int_to_char
!
!
IF ( int < 10 ) THEN
!
IMPLICIT NONE
WRITE( UNIT = int_to_char , FMT = "(I1)" ) int
!
INTEGER, INTENT(IN) :: int
CHARACTER (LEN=6) :: int_to_char
ELSE IF ( int < 100 ) THEN
!
!
IF ( int < 10 ) THEN
!
WRITE( UNIT = int_to_char , FMT = "(I1)" ) int
!
ELSE IF ( int < 100 ) THEN
!
WRITE( UNIT = int_to_char , FMT = "(I2)" ) int
!
ELSE IF ( int < 1000 ) THEN
!
WRITE( UNIT = int_to_char , FMT = "(I3)" ) int
!
ELSE IF ( int < 10000 ) THEN
!
WRITE( UNIT = int_to_char , FMT = "(I4)" ) int
!
ELSE
!
WRITE( UNIT = int_to_char , FMT = "(I5)" ) int
!
END IF
WRITE( UNIT = int_to_char , FMT = "(I2)" ) int
!
RETURN
ELSE IF ( int < 1000 ) THEN
!
END FUNCTION int_to_char
WRITE( UNIT = int_to_char , FMT = "(I3)" ) int
!
ELSE IF ( int < 10000 ) THEN
!
WRITE( UNIT = int_to_char , FMT = "(I4)" ) int
!
ELSE
!
WRITE( UNIT = int_to_char , FMT = "(I5)" ) int
!
END IF
!
RETURN
!
END FUNCTION int_to_char
!
!
!--------------------------------------------------------------------------
SUBROUTINE delete_if_present( filename )
!--------------------------------------------------------------------------
!
IMPLICIT NONE
!
CHARACTER(LEN=*) :: filename
LOGICAL :: exst, opnd
INTEGER :: iunit
!
INQUIRE( FILE = filename, EXIST = exst )
!
IF ( exst ) THEN
!
unit_loop: DO iunit = 99, 1, - 1
!
INQUIRE( UNIT = iunit, OPENED = opnd )
!
IF ( .NOT. opnd ) THEN
!
OPEN( UNIT = iunit, FILE = filename , STATUS = 'OLD' )
CLOSE( UNIT = iunit, STATUS = 'DELETE' )
WRITE( UNIT = stdout, &
FMT = '(/,5X"WARNING: ",A, &
&" file was present; old file deleted")' ) filename
!
RETURN
!
END IF
!
END DO unit_loop
!
CALL errore( 'delete_if_present', 'free unit not found ?!?', 1 )
!
END IF
!
RETURN
!
END SUBROUTINE
!
!
!-----------------------------------------------------------------------
logical function matches (string1, string2)
!-----------------------------------------------------------------------

2
Modules/read_cards.f90
View File

@ -171,7 +171,7 @@ MODULE read_cards_module
ELSE IF ( TRIM(card) == 'CONSTRAINTS' ) THEN
!
CALL card_constraints( input_line )
IF ( ( prog == 'PW' .OR. prog == 'CP' ) .AND. ionode ) &
IF ( ( prog == 'CP' ) .AND. ionode ) &
WRITE( stdout,'(A)') 'Warning: card '//trim(input_line)//' ignored'
!
ELSE IF ( TRIM(card) == 'VHMEAN' ) THEN

2
Modules/read_namelists.f90
View File

@ -282,7 +282,7 @@ MODULE read_namelists_module
!
!
! ... ( 'sd' | 'cg' | 'damp' | 'verlet' | 'none' )
! ... ( 'constrained-verlet' | 'bfgs' | 'constrained-bfgs' | 'beeman' )
! ... ( 'constrained-verlet' | 'bfgs' | 'constrained-damp' | 'beeman' )
!
ion_dynamics = 'none'
ion_radius = 0.5D0

5
PW/Makefile
View File

@ -56,6 +56,7 @@ clean_pw.o \
close_files.o \
compute_dip.o \
constrain.o \
constrains_module.o \
conv_to_num.o \
coset.o \
cryst_to_car.o \
@ -87,7 +88,6 @@ error.o \
error_handler.o \
estimate.o \
ewald.o \
ewald_dipole.o \
fft_scatter.o \
fftw.o \
force_cc.o \
@ -117,8 +117,6 @@ hinit0.o \
hinit1.o \
io_routines.o \
init_ns.o \
init_paw_1.o \
init_paw_2.o \
init_pool.o \
init_run.o \
init_us_1.o \
@ -241,7 +239,6 @@ swap.o \
symrho.o \
symtns.o \
symvect.o \
symz.o \
tabd.o \
trntns.o \
trnvecc.o \

75
PW/add_vuspsi.f90
View File

@ -21,21 +21,28 @@ SUBROUTINE add_vuspsi( lda, n, m, psi, hpsi )
! m number of states psi
! psi
! output:
! hpsi V_US*psi is added to hpsi
! hpsi V_US|psi> is added to hpsi
!
!
USE kinds, ONLY : DP
USE kinds, ONLY : DP
USE basis, ONLY : nat, ntyp, ityp
USE lsda_mod, ONLY : current_spin
USE wvfct, ONLY : gamma_only
USE us, ONLY : vkb, nkb, nh, deeq
!
IMPLICIT NONE
!
!
! ... I/O variables
!
INTEGER, INTENT(IN) :: lda, n, m
COMPLEX(KIND=DP), INTENT(IN) :: psi(lda,m)
COMPLEX(KIND=DP), INTENT(OUT) :: hpsi(lda,m)
!
! ... here the local variables
!
INTEGER :: jkb, ikb, ih, jh, na, nt, ijkb0, ibnd
! counters
!
!
CALL start_clock( 'add_vuspsi' )
!
@ -63,37 +70,48 @@ SUBROUTINE add_vuspsi( lda, n, m, psi, hpsi )
!
IMPLICIT NONE
!
! ... here the local variables
!
INTEGER :: jkb, ikb, ih, jh, na, nt, ijkb0, ibnd
! counters
!
!
IF ( nkb == 0 ) RETURN
!
becp_(:,:) = 0.D0
!
ijkb0 = 0
!
DO nt = 1, ntyp
!
DO na = 1, nat
!
IF ( ityp(na) == nt ) THEN
!
DO ibnd = 1, m
!
DO jh = 1, nh(nt)
!
jkb = ijkb0 + jh
!
DO ih = 1, nh(nt)
!
ikb = ijkb0 + ih
becp_(ikb, ibnd) = becp_(ikb, ibnd) + &
!
becp_(ikb,ibnd) = becp_(ikb,ibnd) + &
deeq(ih,jh,na,current_spin) * becp(jkb,ibnd)
!
END DO
!
END DO
!
END DO
!
ijkb0 = ijkb0 + nh(nt)
!
END IF
!
END DO
!
END DO
!
CALL DGEMM( 'N', 'N', 2*n, m, nkb, 1.d0, vkb, &
2*lda, becp_, nkb, 1.d0, hpsi, 2*lda )
CALL DGEMM( 'N', 'N', ( 2 * n ), m, nkb, 1.D0, vkb, &
( 2 * lda ), becp_, nkb, 1.D0, hpsi, ( 2 * lda ) )
!
RETURN
!
@ -110,38 +128,53 @@ SUBROUTINE add_vuspsi( lda, n, m, psi, hpsi )
!
! ... here the local variables
!
INTEGER :: jkb, ikb, ih, jh, na, nt, ijkb0, ibnd
! counters
COMPLEX(KIND=DP), ALLOCATABLE :: ps (:,:)
! the product vkb and psi
COMPLEX(KIND=DP), ALLOCATABLE :: ps(:,:)
! the product vkb and psi
!
!
IF ( nkb == 0 ) RETURN
!
ALLOCATE( ps( nkb, m) )
!
ALLOCATE( ps( nkb, m ) )
!
ps(:,:) = ( 0.D0, 0.D0 )
!
ijkb0 = 0
!
DO nt = 1, ntyp
!
DO na = 1, nat
!
IF ( ityp(na) == nt ) THEN
!
DO ibnd = 1, m
!
DO jh = 1, nh(nt)
!
jkb = ijkb0 + jh
!
DO ih = 1, nh(nt)
!
ikb = ijkb0 + ih
ps(ikb, ibnd) = ps(ikb, ibnd) + &
!
ps(ikb,ibnd) = ps(ikb,ibnd) + &
deeq(ih,jh,na,current_spin) * becp(jkb,ibnd)
!
END DO
!
END DO
!
END DO
!
ijkb0 = ijkb0 + nh(nt)
!
END IF
!
END DO
!
END DO
!
CALL ZGEMM( 'N', 'N', n, m, nkb, (1.D0, 0.D0) , vkb, &
lda, ps, nkb, (1.D0, 0.D0) , hpsi, lda )
CALL ZGEMM( 'N', 'N', n, m, nkb, ( 1.D0, 0.D0 ) , vkb, &
lda, ps, nkb, ( 1.D0, 0.D0 ) , hpsi, lda )
!
DEALLOCATE( ps )
!
@ -149,4 +182,4 @@ SUBROUTINE add_vuspsi( lda, n, m, psi, hpsi )
!
END SUBROUTINE add_vuspsi_k
!
END SUBROUTINE add_vuspsi
END SUBROUTINE add_vuspsi

6
PW/clean_pw.f90
View File

@ -36,6 +36,7 @@ SUBROUTINE clean_pw()
USE afftnec, ONLY : auxp
#endif
USE fft_types, ONLY : fft_dlay_deallocate
USE constrains_module, ONLY : constr, target
!
IMPLICIT NONE
!
@ -138,6 +139,11 @@ SUBROUTINE clean_pw()
!
CALL berry_closeup( )
!
! ... vectors for ionic constrains
!
IF ( ALLOCATED( constr ) ) DEALLOCATE( constr )
IF ( ALLOCATED( target ) ) DEALLOCATE( target )
!
RETURN
!
END SUBROUTINE clean_pw

877
PW/dynamics.f90
View File

@ -1,468 +1,465 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2004 PWSCF group
! 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 .
!
!-----------------------------------------------------------------------
subroutine dynamics
!-----------------------------------------------------------------------
!
! This routine performs one step of molecular dynamics evolution using
! the Verlet algorithm. Parameters:
! mass mass of the atoms
! dt time step
! temperature starting temperature
! The starting velocities of atoms are set accordingly
! to the starting temperature, in random directions.
! The initial velocity distribution is therefore a constant
! delta_T, nraise are used to change the temperature as follows:
! delta_T = 1 : nothing is done.
! delta_T <> 1 and delta_T >0 : every step the actual temperature is
! multiplied by delta_T, this is done
! rescaling all the velocities
! delta_T < 0 : every 'nraise' step the temperature
! reduced by -delta_T
!
! DA 1997
!
#include "machine.h"
USE io_global, ONLY : stdout
USE kinds, ONLY: DP
USE constants, ONLY: amconv
USE basis, ONLY: nat, ntyp, tau, ityp, atm
USE brilz, ONLY: alat
USE dynam, ONLY: amass, temperature, dt, delta_t, nraise
USE ener, ONLY: etot
USE force_mod, ONLY: force
USE klist, ONLY: nelec
USE relax, ONLY: fixatom, dtau_ref, starting_diag_threshold
USE control_flags, ONLY: ethr, upscale, tr2, imix, alpha0, beta0, istep
use io_files, only : prefix
#ifdef __PARA
use para
!
!----------------------------------------------------------------------------
SUBROUTINE dynamics()
!----------------------------------------------------------------------------
!
! ... This routine performs one step of molecular dynamics evolution using
! ... the Velocity Verlet algorithm.
!
! ... Parameters:
! ... mass mass of the atoms
! ... dt time step
! ... temperature starting temperature
! ... The starting velocities of atoms are set accordingly
! ... to the starting temperature, in random directions.
! ... The initial velocity distribution is therefore a constant
!
! ... delta_T, nraise are used to change the temperature as follows:
!
! ... delta_T = 1 : nothing is done.
! ... delta_T <> 1 and delta_T > 0 : every step the actual temperature is
! ... multiplied by delta_T, this is done
! ... rescaling all the velocities
! ... delta_T < 0 : every 'nraise' step the temperature
! ... reduced by -delta_T
!
! ... Dario Alfe 1997 and Carlo Sbraccia 2004
!
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
USE constants, ONLY : amconv, eps8
USE basis, ONLY : nat, ntyp, tau, ityp, atm
USE brilz, ONLY : alat
USE dynam, ONLY : amass, temperature, dt, delta_t, nraise
USE ener, ONLY : etot
USE force_mod, ONLY : force
USE klist, ONLY : nelec
USE relax, ONLY : if_pos, fixatom, epse, epsf
USE control_flags, ONLY : ethr, upscale, tr2, imix, alpha0, beta0, istep, &
lconstrain, ldamped, conv_ions
USE io_files, ONLY : prefix
#if defined (__PARA)
USE para, ONLY : me, mypool
#endif
implicit none
integer :: natoms ! number of moving atoms
real(kind=DP), allocatable :: &
tauold (:,:,:) ! previous positions of atoms
real(kind=DP), allocatable :: &
a (:,:), mass (:) ! accelerations and masses of atoms
real(kind=DP) :: ekin ! ionic kinetic energy
real(kind=DP) :: total_mass, temp_new, tempo, norm_of_dtau
real(kind=DP) :: ml (3) ! total linear momentum
integer :: na, ipol, it ! counters
logical :: exst
real(kind=DP), parameter :: convert_E_to_temp= 315642.28d0 * 0.5d0, &
eps = 1.d-6
allocate (mass( nat))
allocate (a( 3, nat))
allocate (tauold( 3, nat, 3))
tauold(:,:,:) = 0.d0
dtau_ref = 0.2
natoms = nat - fixatom
!
! one Ryd a.u. of time is 4.84*10^-17 seconds, i.e. 0.0484 femtosecond
IMPLICIT NONE
!
call seqopn (4, trim(prefix)//'.md', 'formatted', exst)
if (.not.exst) then
close (unit = 4, status = 'delete')
WRITE( stdout, '(/5x,"Starting temperature = ",f8.2," K")') &
temperature
do na = 1, ntyp
WRITE( stdout, '(5x,"amass(",i1,") = ",f6.2)') na, amass (na)
enddo
WRITE( stdout, '(5x,"Time step = ",f6.2," a.u., ",f6.4, &
& " femto-seconds")') dt, dt * 0.0484
! ... local variables
!
REAL(KIND=DP), ALLOCATABLE :: vel(:,:), acc(:,:), accold(:,:)
! velocities, accelerations and accelerations of the previous step
REAL(KIND=DP), ALLOCATABLE :: mass(:)
! masses of atoms
REAL(KIND=DP) :: ekin, etotold
! ionic kinetic energy
! ionic potential energy (of the previous step)
REAL(KIND=DP) :: total_mass, temp_new, elapsed_time, norm_of_dtau
REAL(KIND=DP) :: ml(3), mlt
! total linear momentum and its modulus
INTEGER :: i, na, it
! counters
LOGICAL :: exst
REAL(KIND=DP), PARAMETER :: convert_E_to_temp = 315642.28D0 * 0.5D0
!
!
ALLOCATE( mass( nat ) )
ALLOCATE( vel( 3, nat ) )
ALLOCATE( acc( 3, nat ) )
ALLOCATE( accold( 3, nat ) )
!
vel = 0.D0
acc = 0.D0
accold = 0.D0
!
! ... one Ryd a.u. of time is 4.84*10^-17 seconds, i.e. 0.0484 femtoseconds
!
CALL seqopn( 4, TRIM( prefix ) // '.md', 'FORMATTED', exst )
!
IF ( .NOT. exst ) THEN
!
! masses in atomic rydberg units
CLOSE( UNIT = 4, STATUS = 'DELETE' )
!
total_mass = 0.d0
do na = 1, nat
mass (na) = amass (ityp (na) ) * amconv
total_mass = total_mass + mass (na)
enddo
WRITE( stdout, '(/,5X,"Starting temperature = ",F8.2," K")' ) &
temperature
!
! initial thermalization. N.B. tau is in units of alat
DO na = 1, ntyp
!
WRITE( stdout, '(5X,"amass(",I1,") = ",F6.2)' ) na, amass(na)
!
END DO
!
WRITE( stdout, '(5X,"Time step = ",F6.2," a.u., ",F6.4, &
& " femto-seconds")' ) dt, ( dt * 0.0484D0 )
!
!... masses in atomic rydberg units
!
total_mass = 0.D0
!
DO na = 1, nat
!
mass(na) = amass( ityp(na) ) * amconv
total_mass = total_mass + mass(na)
!
END DO
!
! ... initial thermalization. N.B. tau is in units of alat
!
CALL start_therm()
!
elapsed_time = 0.D0
!
call start_therm (mass, tauold)
tempo = 0.d0
temp_new = temperature
!
it = 0
else
read (4, * ) temp_new, mass, total_mass, tauold, tempo, it
close (unit = 4, status = 'keep')
!
ELSE
!
READ( UNIT = 4, FMT = * ) &
etotold, temp_new, mass, total_mass, elapsed_time, it, tau, vel, accold
CLOSE( UNIT = 4, STATUS = 'KEEP' )
!
istep = it + 1
endif
tempo = tempo + dt * 0.0000484
!
END IF
!
elapsed_time = elapsed_time + dt * 0.0000484D0
!
it = it + 1
if (mod (it, nraise) .eq.0.and.delta_T.lt.0) then
WRITE( stdout, '(/5x,"Thermalization: delta_T = ",f6.3, &
& ", T = ",f6.1)') - delta_T, temp_new - delta_T
call thermalize (temp_new, temp_new - delta_T, tauold)
endif
if (delta_T.ne.1.d0.and.delta_T.ge.0) then
WRITE( stdout, '(/5x,"Thermalization: delta_T = ",f6.3, &
& ", T = ",f6.1)') delta_T, temp_new * delta_T
call thermalize (temp_new, temp_new * delta_T, tauold)
endif
WRITE( stdout, '(/5x,"Entering Dynamics; it = ",i5," time = ", &
& f8.5," pico-seconds"/)') it, tempo
!
! calculate accelerations in a.u. units / alat
IF ( MOD( it, nraise ) == 0 .AND. delta_T < 0 ) THEN
!
WRITE( stdout, '(/,5X,"Thermalization: delta_T = ",F6.3, &
& ", T = ",F6.1)' ) - delta_T, ( temp_new - delta_T )
!
CALL thermalize( temp_new, ( temp_new - delta_T ) )
!
ELSE IF ( delta_T /= 1.D0 .AND. delta_T >= 0 ) THEN
!
WRITE( stdout, '(/,5X,"Thermalization: delta_T = ",F6.3, &
& ", T = ",F6.1)' ) delta_T, temp_new * delta_T
!
CALL thermalize( temp_new, temp_new * delta_T )
!
END IF
!
do na = 1, natoms
a (1, na) = force (1, na) / mass (na) / alat
a (2, na) = force (2, na) / mass (na) / alat
a (3, na) = force (3, na) / mass (na) / alat
enddo
! ... check if convergence for structural minimization is achieved
!
! save the previous two steps ( a total of three steps is saved)
!
tauold (:, :, 3) = tauold (:, :, 2)
tauold (:, :, 2) = tauold (:, :, 1)
!
! move the atoms accordingly to the classical equation of motion
!
call verlet (tau, tauold, a, natoms, ekin, mass, ml, dt)
!
! find the best coefficients for the extrapolation of the potential
!
call find_alpha_and_beta (nat, tau, tauold, alpha0, beta0)
#ifdef __PARA
if (me.eq.1) call poolbcast (1, alpha0)
call broadcast (1, alpha0)
if (me.eq.1) call poolbcast (1, beta0)
call broadcast (1, beta0)
#endif
!
! calculate the "norm" of the step (used to update the convergence threshold)
!
norm_of_dtau = 0.d0
do na = 1, nat
do ipol = 1, 3
norm_of_dtau = norm_of_dtau + (tau(ipol,na) - tauold(ipol,na,1)) **2
enddo
enddo
norm_of_dtau = sqrt (norm_of_dtau)
!
! ... ethr is now updated in electrons with a different procedure
! ... this value of ethr is overwritten apart when the old style
! ... update is used (see OLDSTYLE precompiler variable in electrons)
!
if (imix.lt.0) then
ethr = starting_diag_threshold * &
max (1.d0 / upscale, min (1.d0, norm_of_dtau / dtau_ref) ) **2
else
ethr = tr2 / nelec
end if
!
! find the new temperature
!
temp_new = 2.d0 / 3.d0 * ekin * alat**2 / natoms * convert_E_to_temp
!
! save on file needed quantity
!
call seqopn (4, trim(prefix)//'.md', 'formatted', exst)
write (4, * ) temp_new, mass, total_mass, tauold, tempo, it
close (unit = 4, status = 'keep')
!
call output_tau( .FALSE. )
WRITE( stdout, '(/5x,"Ekin = ",f14.8," Ryd T = ",f6.1," K ", &
& " Etot = ",f14.8)') ekin*alat**2, temp_new, ekin*alat**2+etot
!
! total linear momentum must be zero if all atoms move
!
if (fixatom == 0) then
if ( abs(ml(1)) > eps .OR. abs(ml(2)) > eps .OR. abs(ml(3)) > eps ) then
call errore ('dynamics', 'Total linear momentum <> 0', - 1)
WRITE( stdout, '(5x,"Linear momentum: ",3f12.8)') ml
end if
endif
deallocate (tauold)
deallocate (a)
deallocate (mass)
return
end subroutine dynamics
!---------------------------------------------------------------------
subroutine verlet (rnew, rold, a, n, ec, mass, ml, dt)
!---------------------------------------------------------------------
!
! Verlet algorithm to update atomic position
!
USE kinds, ONLY: DP
implicit none
! INPUT
integer :: n ! number of particles
real(kind=DP) :: dt, & ! time step
a (3, n), & ! accelerations
mass (n) ! atom masses
! OUTPUT
real(kind=DP) :: ec, ml (3) ! kinetic energy and total linear momentum
! INPUT/OUTPUT
real(kind=DP) :: rold (3, n), & ! in: previous, out: present atomic positions
rnew (3, n) ! in: present, out: new atomic positions
! LOCAL
integer :: i
real(kind=DP) :: dtsq, dt2, r (3), v (3)
! dtsq=dt**2, dt2=2*dt
!
ml(:) = 0.d0
ec = 0.d0
dtsq = dt**2
!
dt2 = dt * 2.0
do i = 1, n
r(:) = 2.0 * rnew(:,i) - rold(:,i) + a(:,i) * dtsq
v(:) = (r(:) - rold(:,i) ) / dt2
rold(:,i) = rnew(:,i)
rnew(:,i) = r(:)
ml(:) = ml(:) + v(:) * mass(i)
ec = ec + 0.5 * mass(i) * (v(1)**2 + v(2)**2 + v(3)**2)
enddo
return
end subroutine verlet
!-----------------------------------------------------------------------
subroutine start_therm (mass, tauold)
!-----------------------------------------------------------------------
!
! Starting thermalization of the system
!
USE kinds, ONLY: DP
USE basis, ONLY: nat, tau
USE brilz, ONLY: alat
USE dynam, ONLY: temperature, dt!, delta_t, nraise
USE relax, ONLY: fixatom
USE symme, only: invsym, nsym, irt
#ifdef __PARA
use para
#endif
implicit none
real(kind=DP) :: mass (nat), tauold (3, nat)
!
integer :: na, nb, natoms
real(kind=DP) :: total_mass, temp_new, aux, convert_E_to_temp, velox,&
ek, ml(3), direzione_x, direzione_y, direzione_z, modulo, rndm
! ek = kinetic energy
real(kind=DP), allocatable :: step(:,:)
external rndm
parameter (convert_E_to_temp = 315642.28d0 * 0.5d0)
!
allocate (step(3,nat))
aux = temperature / convert_E_to_temp
natoms = nat - fixatom
ml(:) = 0.d0
!
! velocity in random direction, with modulus accordingly to mass and
! temperature: 3/2KT = 1/2mv^2
!
#ifdef __PARA
!
! only the first processor calculates ...
!
if (me.eq.1.and.mypool.eq.1) then
#endif
do na = 1, natoms
IF ( ldamped ) THEN
!
conv_ions = ( etotold - etot ) < epse
!
DO i = 1, 3
DO na = 1, nat
conv_ions = conv_ions .AND. ( ABS( force(i,na) ) < epsf )
END DO
END DO
!
IF ( conv_ions ) THEN
!
! N.B. velox is in a.u. units /alat
WRITE( UNIT = stdout, &
FMT = '(/,5X,"Damped Dynamics: convergence achieved in ",I3, &
& " steps")' ) it
WRITE( UNIT = stdout, &
FMT = '(/,5X,"Efinal = ",F15.8,/)' ) etot
!
velox = sqrt (3.d0 * aux / mass (na) ) / alat
direzione_x = rndm () - .5d0
direzione_y = rndm () - .5d0
direzione_z = rndm () - .5d0
modulo = sqrt (direzione_x**2 + direzione_y**2 + direzione_z**2)
step (1, na) = velox / modulo * direzione_x
step (2, na) = velox / modulo * direzione_y
step (3, na) = velox / modulo * direzione_z
enddo
#ifdef __PARA
CALL output_tau( .TRUE. )
!
RETURN
!
END IF
!
! ... and distributes the velocities
END IF
!
WRITE( stdout, '(/,5X,"Entering Dynamics; it = ",I5," time = ", &
&F8.5," pico-seconds",/)' ) it, elapsed_time
!
! ... calculate accelerations in a.u. units / alat
!
FORALL( na = 1 : nat ) &
acc(:,na) = force(:,na) / mass(na) / alat
!
! ... atoms are moved accordingly to the classical equation of motion.
! ... Damped dynamics ( based on the quick-min algorithm ) is also
! ... done here.
!
vel = vel + 0.5D0 * dt * ( accold + acc )
!
IF ( ldamped ) CALL projec_velocity()
!
! ... constrains ( atoms kept fixed ) are reinforced
!
vel =vel * REAL( if_pos )
!
tau = tau + dt * vel + 0.5D0 * dt**2 * acc
!
ml = 0.D0
ekin = 0.D0
!
DO na = 1, nat
!
ml(:) = ml(:) + vel(:,na) * mass(na)
ekin = ekin + &
0.5D0 * mass(na) * ( vel(1,na)**2 + vel(2,na)**2 + vel(3,na)**2 )
!
endif
if (me.eq.1) call poolbcast (3 * natoms, step)
call broadcast (3 * natoms, step)
END DO
!
! ... find the new temperature
!
temp_new = 2.D0 / 3.D0 * ekin * alat**2 / nat * convert_E_to_temp
!
! ... save on file needed quantity
!
CALL seqopn( 4, TRIM( prefix ) // '.md', 'FORMATTED', exst )
!
WRITE( UNIT = 4, FMT = * ) &
etot, temp_new, mass, total_mass, elapsed_time, it, tau, vel, acc
!
CLOSE( UNIT = 4, STATUS = 'KEEP' )
!
DO na = 1, nat
!
WRITE( stdout, '(A3,3F12.7)') atm(ityp(na)), tau(:,na)
!
END DO
!
IF ( it == 1 ) THEN
!
WRITE( stdout, '(/,5X,"Ekin = ",F14.8," Ryd T = ",F6.1," K ", &
& " Etot = ",F14.8)' ) &
temperature * 3.D0 / 2.D0 * nat / convert_E_to_temp, temperature, &
( temperature * 3.D0 / 2.D0 * nat / convert_E_to_temp + etot )
!
ELSE
!
WRITE( stdout, '(/,5X,"Ekin = ",F14.8," Ryd T = ",F6.1," K ", &
& " Etot = ",F14.8)' ) &
ekin*alat**2, temp_new, ( ekin*alat**2 + etot )
!
END IF
!
! ... total linear momentum must be zero if all atoms move
!
mlt = ABS( ml(1) ) + ABS( ml(2) ) + ABS( ml(3) )
!
IF ( .NOT. ( ldamped .OR. lconstrain ) .AND. &
( fixatom == 0 ) .AND. ( mlt > eps8 ) ) &
CALL errore( 'dynamics', 'Total linear momentum <> 0', - 1 )
!
WRITE( stdout, '(5X,"Linear momentum: ",3F18.14)' ) ml
!
DEALLOCATE( mass )
DEALLOCATE( vel )
DEALLOCATE( acc )
DEALLOCATE( accold )
!
RETURN
!
CONTAINS
!
! ... internal procedure
!
!-----------------------------------------------------------------------
SUBROUTINE projec_velocity()
!-----------------------------------------------------------------------
!
USE constants, ONLY : eps32
USE basic_algebra_routines
!
IMPLICIT NONE
!
! ... local variables
!
REAL (KIND=DP) :: norm_acc, acc_versor(3,nat)
!
! ... external functions
!
REAL (KIND=DP), EXTERNAL :: DNRM2, DDOT
!
norm_acc = DNRM2( 3*nat, acc, 1 )
!
IF ( norm_acc >= eps32 ) THEN
!
acc_versor = acc / norm_acc
!
vel = MAX( 0.D0, DDOT( 3*nat, vel, 1, acc_versor, 1 ) ) * acc_versor
!
ELSE
!
vel = 0.D0
!
END IF
!
END SUBROUTINE projec_velocity
!
!
!-----------------------------------------------------------------------
SUBROUTINE start_therm()
!-----------------------------------------------------------------------
!
! ... Starting thermalization of the system
!
USE symme, ONLY : invsym, nsym, irt
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: na, nb
REAL(KIND=DP) :: total_mass, aux, velox, ek, &
ml(3), dir_x, dir_y, dir_z, module
! ek = kinetic energy
REAL(KIND=DP),EXTERNAL :: rndm
!
!
aux = temperature / convert_E_to_temp
!
! ... velocity in random direction, with modulus accordingly to mass
! ... and temperature: 3/2KT = 1/2mv^2
!
#if defined (__PARA)
!
! ... only the first processor calculates ...
!
IF ( me == 1 .AND. mypool == 1 ) THEN
#endif
!
! if there is inversion symmetry equivalent atoms have opposite velocities
!
if (invsym) then
do na = 1, natoms
nb = irt (nsym / 2 + 1, na)
if (nb.gt.na) then
step (1, nb) = - step (1, na)
step (2, nb) = - step (2, na)
step (3, nb) = - step (3, na)
endif
!
! the atom on the inversion center is kept fixed
!
if (na.eq.nb) then
step (1, na) = 0.d0
step (2, na) = 0.d0
step (3, na) = 0.d0
endif
enddo
else
!
DO na = 1, nat
!
! ... N.B. velox is in a.u. units /alat
!
velox = SQRT( 3.D0 * aux / mass(na) ) / alat
!
dir_x = rndm() - 0.5D0
dir_y = rndm() - 0.5D0
dir_z = rndm() - 0.5D0
!
module = SQRT( dir_x**2 + dir_y**2 + dir_z**2 )
!
vel(1,na) = velox * dir_x / module
vel(2,na) = velox * dir_y / module
vel(3,na) = velox * dir_z / module
!
END DO
#ifdef __PARA
!
! ... and distributes the velocities
!
END IF
!
IF ( me == 1 ) CALL poolbcast( 3 * nat, vel )
!
CALL broadcast( 3 * nat, vel )
#endif
!
! ... if there is inversion symmetry equivalent atoms have
! ... opposite velocities
!
ml(:) = 0.D0
!
IF ( invsym ) THEN
!
DO na = 1, nat
!
nb = irt( nsym / 2 + 1, na )
!
IF ( nb > na ) vel(:,nb) = - vel(:,na)
!
! ... the atom on the inversion center is kept fixed
!
IF ( na == nb ) vel(:,na) = 0.D0
!
END DO
!
ELSE
!
! ... put total linear momentum equal zero if all atoms move
!
PRINT *, "FIXATOM = ", FIXATOM
IF ( fixatom == 0 ) THEN
!
total_mass = 0.D0
!
DO na = 1, nat
!
total_mass = total_mass + mass(na)
!
ml(:) = ml(:) + mass(na) * vel(:,na)
!
END DO
!
ml(:) = ml(:) / total_mass
!
END IF
!
END IF
!
ek = 0.D0
!
DO na = 1, nat
!
vel(:,na) = vel(:,na) - ml(:)
!
ek = ek + 0.5D0 * mass(na) * ( ( vel(1,na) - ml(1) )**2 + &
( vel(2,na) - ml(2) )**2 + &
( vel(3,na) - ml(3) )**2 )
!
END DO
!
! ... after the velocity of the center of mass has been subtracted the
! ... temperature is usually changed. Set again the temperature to the
! ... right value.
!
temp_new = 2.D0 * ek / ( 3.D0 * nat ) * alat**2 * convert_E_to_temp
!
CALL thermalize( temp_new, temperature )
!
RETURN
!
END SUBROUTINE start_therm
!
! put total linear momentum equal zero if all atoms move
!
if (fixatom.eq.0) then
total_mass = 0.d0
do na = 1, natoms
total_mass = total_mass + mass (na)
ml (1) = ml (1) + mass (na) * step (1, na)
ml (2) = ml (2) + mass (na) * step (2, na)
ml (3) = ml (3) + mass (na) * step (3, na)
enddo
ml (1) = ml (1) / total_mass
ml (2) = ml (2) / total_mass
ml (3) = ml (3) / total_mass
endif
endif
!
! -step is the velocity
!
ek = 0.d0
do na = 1, natoms
tauold (1, na) = (step (1, na) - ml (1) ) * dt + tau (1, na)
tauold (2, na) = (step (2, na) - ml (2) ) * dt + tau (2, na)
tauold (3, na) = (step (3, na) - ml (3) ) * dt + tau (3, na)
ek = ek + 0.5d0 * mass (na) * ( (step (1, na) - ml (1) ) **2 + &
(step (2, na) - ml (2) ) **2 + (step (3, na) - ml (3) ) **2)
enddo
!
! after the velocity of the center of mass has been subtracted the
! temperature is usually changed. Set again the temperature to the
! right value.
!
temp_new = 2.d0 * ek / (3.d0 * natoms) * alat**2 * convert_E_to_temp
call thermalize (temp_new, temperature, tauold)
deallocate(step)
return
end subroutine start_therm
!-------------------------------------------------------------------
subroutine thermalize (temp_old, temp_new, tauold)
!-------------------------------------------------------------------
USE kinds, ONLY: DP
USE basis, ONLY: nat, tau
USE relax, ONLY: fixatom
USE dynam, ONLY: dt
implicit none
real(kind=DP) :: tauold (3, nat), temp_new, temp_old
!
integer :: na, natoms
real(kind=DP) :: velox, aux
!
natoms = nat - fixatom
!
! rescale the velocities by a factor 3/2KT/Ek
!
if (temp_new.gt.0.d0.and.temp_old.gt.0.d0) then
aux = sqrt (temp_new / temp_old)
else
aux = 0.d0
endif
do na = 1, natoms
velox = (tau (1, na) - tauold (1, na) ) / dt
tauold (1, na) = tau (1, na) - dt * aux * velox
velox = (tau (2, na) - tauold (2, na) ) / dt
tauold (2, na) = tau (2, na) - dt * aux * velox
velox = (tau (3, na) - tauold (3, na) ) / dt
tauold (3, na) = tau (3, na) - dt * aux * velox
enddo
!
return
end subroutine thermalize
!-----------------------------------------------------------------------
subroutine find_alpha_and_beta (nat, tau, tauold, alpha0, beta0)
!-----------------------------------------------------------------------
!
! This routine finds the best coefficients alpha0 and beta0 so that
!
! | tau(t+dt) - tau' | is minimum, where
!
! tau' = alpha0 * ( tau(t) - tau(t-dt) ) +
! beta0 * ( tau(t-dt) - tau(t-2*dt) )
!
USE kinds, ONLY: DP
implicit none
integer :: nat, na, ipol
real(kind=DP) :: chi, alpha0, beta0, tau (3, nat), tauold (3, nat, 3)
real(kind=DP) :: a11, a12, a21, a22, b1, b2, c, det
!
! solution of the linear system
!
a11 = 0.d0
a12 = 0.d0
a21 = 0.d0
a22 = 0.d0 + 1.d-12
b1 = 0.d0
b2 = 0.d0
c = 0.d0
do na = 1, nat
do ipol = 1, 3
a11 = a11 + (tauold (ipol, na, 1) - tauold (ipol, na, 2) ) **2
a12 = a12 + (tauold (ipol, na, 1) - tauold (ipol, na, 2) ) &
* (tauold (ipol, na, 2) - tauold (ipol, na, 3) )
a22 = a22 + (tauold (ipol, na, 2) - tauold (ipol, na, 3) ) **2
b1 = b1 - (tauold (ipol, na, 1) - tau (ipol, na) ) * (tauold ( &
ipol, na, 1) - tauold (ipol, na, 2) )
b2 = b2 - (tauold (ipol, na, 1) - tau (ipol, na) ) * (tauold ( &
ipol, na, 2) - tauold (ipol, na, 3) )
c = c + (tauold (ipol, na, 1) - tau (ipol, na) ) **2
enddo
enddo
a21 = a12
!
det = a11 * a22 - a12 * a21
if (det.lt.0d0) call errore ('find_alpha_and_beta', ' det.le.0', 1)
!
! det > 0 case: a well defined minimum exists
!
if (det.gt.0d0) then
alpha0 = (b1 * a22 - b2 * a12) / det
beta0 = (a11 * b2 - a21 * b1) / det
else
!
! det = 0 case: the two increments are linearly dependent, chose
! solution with beta=0 (discard oldest configuration)
!
alpha0 = 1.d0
if (a11.gt.0.d0) alpha0 = b1 / a11
beta0 = 0.d0
endif
chi = 0.d0
do na = 1, nat
do ipol = 1, 3
chi = chi + ( (1 + alpha0) * tauold (ipol, na, 1) + (beta0 - &
alpha0) * tauold (ipol, na, 2) - beta0 * tauold (ipol, na, 3) &
- tau (ipol, na) ) **2
enddo
enddo
! WRITE( stdout, * ) chi, alpha0, beta0
return
end subroutine find_alpha_and_beta
!-----------------------------------------------------------------------
SUBROUTINE thermalize( temp_old, temp_new )
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!
! ... INPUT variables
!
REAL(KIND=DP) :: temp_new, temp_old
!
! ... local variables
!
INTEGER :: na
REAL(KIND=DP) :: aux
!
!
! ... rescale the velocities by a factor 3 / 2KT / Ek
!
IF ( temp_new > 0.D0 .AND. temp_old > 0.D0 ) THEN
!
aux = SQRT( temp_new / temp_old )
!
ELSE
!
aux = 0.D0
!
END IF
!
vel = vel * aux * REAL( if_pos )
!
RETURN
!
END SUBROUTINE thermalize
!
END SUBROUTINE dynamics

157
PW/input.f90
View File

@ -91,7 +91,8 @@ SUBROUTINE iosys()
nosym_ => nosym, &
modenum_ => modenum, &
reduce_io, ethr, lscf, lbfgs, lmd, lneb, lphonon, &
noinv, time_max, restart, loldbfgs
noinv, time_max, restart, loldbfgs, lconstrain, &
ldamped
USE wvfct, ONLY : ibm_baco2, &
nbnd_ => nbnd
USE fixed_occ, ONLY : tfixed_occ
@ -499,16 +500,23 @@ SUBROUTINE iosys()
& ' fixed occupations and lsda not implemented ', 1 )
END IF
!
calc = ' '
! ... initialization
!
calc = ' '
lbfgs = .FALSE.
loldbfgs = .FALSE.
lmd = .FALSE.
ldamped = .FALSE.
lconstrain = .FALSE.
!
IF ( TRIM( calculation ) == 'relax' ) THEN
!
SELECT CASE ( TRIM( ion_dynamics ) )
CASE ( 'bfgs' )
iswitch = 1 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
lbfgs = .TRUE.
loldbfgs = .FALSE.
epse = etot_conv_thr
epsf = forc_conv_thr
epse = etot_conv_thr
epsf = forc_conv_thr
!
IF ( epse <= 20.D0 * ( tr2 / upscale ) ) &
CALL errore( ' iosys ', ' required etot_conv_thr is too small:' // &
@ -516,25 +524,25 @@ SUBROUTINE iosys()
!
CASE ( 'old-bfgs' )
iswitch = 1 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
lbfgs = .FALSE.
loldbfgs = .TRUE.
epse = etot_conv_thr
epsf = forc_conv_thr
epse = etot_conv_thr
epsf = forc_conv_thr
!
IF ( epse <= 20.D0 * ( tr2 / upscale ) ) &
CALL errore( ' iosys ', ' required etot_conv_thr is too small:' // &
& ' conv_thr must be reduced', 1 )
!
CASE ( 'constrained-bfgs' )
iswitch = 2 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
lbfgs = .TRUE.
epse = etot_conv_thr
epsf = forc_conv_thr
CASE ( 'constrained-damp' )
iswitch = 4 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
lmd = .TRUE.
ldamped = .TRUE.
lconstrain = .TRUE.
epse = etot_conv_thr
epsf = forc_conv_thr
CASE ( 'damp' )
iswitch = 3 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
lbfgs = .FALSE.
lmd = .TRUE.
calc = 'mm'
ldamped = .TRUE.
epse = etot_conv_thr
epsf = forc_conv_thr
ntcheck = nstep + 1
@ -543,13 +551,16 @@ SUBROUTINE iosys()
& ': ion_dynamics=' // TRIM( ion_dynamics ) // &
& ' not supported', 1 )
END SELECT
END IF
!
IF ( TRIM( calculation ) == 'md' ) THEN
!
ELSE IF ( TRIM( calculation ) == 'md' ) THEN
!
lmd = .TRUE.
!
SELECT CASE ( TRIM( ion_dynamics ) )
CASE ( 'verlet' )
iswitch = 3 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
CASE ( 'constrained-verlet' )
lconstrain = .TRUE.
iswitch = 4 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
CASE ( 'beeman' )
iswitch = 3 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
@ -560,9 +571,9 @@ SUBROUTINE iosys()
& ': ion_dynamics=' // TRIM( ion_dynamics ) // &
& ' not supported', 1 )
END SELECT
END IF
!
IF ( TRIM( calculation ) == 'vc-relax' ) THEN
!
ELSE IF ( TRIM( calculation ) == 'vc-relax' ) THEN
!
SELECT CASE ( TRIM( cell_dynamics ) )
CASE ( 'none' )
epse = etot_conv_thr
@ -592,9 +603,9 @@ SUBROUTINE iosys()
& ': ion_dynamics=' // TRIM( ion_dynamics ) // &
& ' not supported', 1 )
END IF
END IF
!
IF ( TRIM( calculation ) == 'vc-md' ) THEN
!
ELSE IF ( TRIM( calculation ) == 'vc-md' ) THEN
!
SELECT CASE ( TRIM( cell_dynamics ) )
CASE ( 'none' )
iswitch = 3 ! ... obsolescent: do not use in new code ( 29/10/2003 C.S.)
@ -613,16 +624,16 @@ SUBROUTINE iosys()
& ': ion_dynamics=' // TRIM( ion_dynamics ) // &
& ' not supported', 1 )
END SELECT
!
IF ( TRIM( ion_dynamics ) /= 'beeman' ) THEN
CALL errore( ' iosys ', 'calculation=' // TRIM( calculation ) // &
& ': ion_dynamics=' // TRIM( ion_dynamics ) // &
& ' not supported', 1 )
END IF
END IF
!
! ... NEB specific
!
IF ( TRIM( calculation ) == 'neb' ) THEN
!
ELSE IF ( TRIM( calculation ) == 'neb' ) THEN
!
! ... NEB specific
!
IF ( num_of_images < 2 ) THEN
CALL errore( ' iosys ', 'calculation=' // TRIM( calculation ) // &
@ -646,38 +657,28 @@ SUBROUTINE iosys()
!
END IF
!
SELECT CASE ( minimization_scheme )
! ... initialization
!
lsteep_des = .FALSE.
lquick_min = .FALSE.
ldamped_dyn = .FALSE.
lmol_dyn = .FALSE.
!
SELECT CASE ( minimization_scheme )
CASE ( "sd" )
lsteep_des = .TRUE.
lquick_min = .FALSE.
ldamped_dyn = .FALSE.
lmol_dyn = .FALSE.
CASE ( "quick-min" )
lsteep_des = .FALSE.
lquick_min = .TRUE.
ldamped_dyn = .FALSE.
lmol_dyn = .FALSE.
CASE ( "damped-dyn" )
lsteep_des = .FALSE.
lquick_min = .FALSE.
ldamped_dyn = .TRUE.
lmol_dyn = .FALSE.
CASE ( "mol-dyn" )
lsteep_des = .FALSE.
lquick_min = .FALSE.
ldamped_dyn = .FALSE.
lmol_dyn = .TRUE.
IF ( temp_req == 0 ) &
WRITE( stdout,'(/,T2,"WARNING: tepm_req has not been set" )')
!
temp_req = temp_req / ( eV_to_kelvin * AU )
!
CASE default
!
CALL errore( ' iosys ','calculation=' // TRIM( calculation ) // &
& ': unknown minimization_scheme', 1 )
!
END SELECT
!
END IF
@ -1030,9 +1031,8 @@ END SUBROUTINE iosys
SUBROUTINE read_cards( psfile, atomic_positions_ )
!-----------------------------------------------------------------------
!
USE parser
USE wvfct, ONLY : gamma_only
USE brilz, ONLY : at, ibrav, symm_type
USE brilz, ONLY : at, ibrav, symm_type, celldm
USE basis, ONLY : nat, ntyp, ityp, tau, atm
USE klist, ONLY : nks
USE ktetra, ONLY : nk1_ => nk1, &
@ -1049,14 +1049,20 @@ SUBROUTINE read_cards( psfile, atomic_positions_ )
USE relax, ONLY : fixatom, &
if_pos_ => if_pos
USE dynam, ONLY : amass
USE control_flags, ONLY : lconstrain
USE constrains_module, ONLY : nconstr, constr_tol, constr, target
USE input_parameters, ONLY : atom_label, atom_pfile, atom_mass, &
atom_ptyp, taspc, tapos, rd_pos, &
atomic_positions, if_pos, sp_pos, &
k_points, xk, wk, nk1, nk2, nk3, &
k1, k2, k3, nkstot, cell_symmetry, rd_ht, &
trd_ht, f_inp, calculation
trd_ht, f_inp, calculation,&
nconstr_inp, constr_tol_inp, constr_inp
USE read_cards_module, ONLY : read_cards_base => read_cards
!
USE parser
USE basic_algebra_routines, ONLY : norm
!
IMPLICIT NONE
!
CHARACTER (LEN=80) :: psfile(ntyp)
@ -1203,7 +1209,36 @@ SUBROUTINE read_cards( psfile, atomic_positions_ )
IF ( ibrav /= 0 .AND. tcell ) &
CALL errore( ' cards ', ' redundant data for cell parameters', 2 )
!
RETURN
IF ( lconstrain ) THEN
!
nconstr = nconstr_inp
constr_tol = constr_tol_inp
!
ALLOCATE( constr(2,nconstr) )
ALLOCATE( target(nconstr) )
!
constr(:,:) = constr_inp(:,1:nconstr)
!
! ... target value of the constrain ( in bohr )
!
DO ia = 1, nconstr
!
target(ia) = norm( tau(:,constr(1,ia)) - &
tau(:,constr(2,ia)) ) * celldm(1)
!
END DO
!
! ... "if_pos" constrains are not compatible with other constrains
!
if_pos_ = 1
!
IF ( fixatom > 0 ) &
CALL errore( ' cards ', &
& ' constrains are not compatible with fixed atoms', 1 )
!
END IF
!
RETURN
!
END SUBROUTINE read_cards
!
@ -1221,7 +1256,7 @@ SUBROUTINE verify_tmpdir()
USE mp, ONLY : mp_barrier
#endif
!
USE parser, ONLY : int_to_char
USE parser, ONLY : int_to_char, delete_if_present
!
IMPLICIT NONE
!
@ -1258,13 +1293,19 @@ SUBROUTINE verify_tmpdir()
!
#if defined (__PARA)
IF ( me == 1 .AND. mypool == 1 ) THEN
#endif
#endif
!
! ... wfc-extrapolation file is removed
!
CALL delete_if_present( TRIM( tmp_dir ) // TRIM( prefix ) // '.update' )
!
! ... MD restart file is removed
!
CALL delete_if_present( TRIM( tmp_dir ) // TRIM( prefix ) // '.md' )
!
! ... BFGS rstart file is removed
!
OPEN( UNIT = 4, FILE = TRIM( tmp_dir ) // TRIM( prefix ) // '.bfgs', &
STATUS = 'UNKNOWN' )
CLOSE( UNIT = 4, STATUS = 'DELETE' )
CALL delete_if_present( TRIM( tmp_dir ) // TRIM( prefix ) // '.bfgs' )
!
#if defined (__PARA)
END IF
@ -1323,11 +1364,9 @@ SUBROUTINE verify_tmpdir()
IF ( me == 1 .AND. mypool == 1 ) THEN
#endif
!
! ... standard output of the self consistency is removed
! ... standard output of the self-consistency is removed
!
OPEN( UNIT = 4, FILE = TRIM( tmp_dir ) // 'PW.out', &
STATUS = 'UNKNOWN' )
CLOSE( UNIT = 4, STATUS = 'DELETE' )
CALL delete_if_present( TRIM( tmp_dir ) // 'PW.out' )
!
#if defined (__PARA)
END IF

494
PW/move_ions.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001-2003 PWSCF group
! Copyright (C) 2001-2004 PWSCF group
! 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,
@ -13,7 +13,7 @@ SUBROUTINE move_ions()
!
! ... This routine moves the ions according to the requested scheme:
!
! ... iswitch = 1 bfgs minimizations or conjugate gradient
! ... iswitch = 1 bfgs minimizations
! ... iswitch = 2 constrained bfgs minimization:
! ... the user must supply the routine 'constrain' which
! ... defines the constraint equation and the gradient
@ -28,12 +28,15 @@ SUBROUTINE move_ions()
! ... dgv(i,na) = ---------------.
! ... D tau(i,na)
!
! ... iswitch = 3 molecular dynamics, (verlet of vcsmd)
! ... iswitch = 3 molecular dynamics, ( verlet of vcsmd )
! ... iswitch = 4 molecular dynamics with one constraint,
! ... the same conventions as iswitch = 2
!
! ... coefficients for potential and wavefunctions extrapolation are
! ... also computed here
!
USE io_global, ONLY : stdout
USE io_files, ONLY : tmp_dir
USE io_files, ONLY : tmp_dir, prefix
USE bfgs_module, ONLY : lbfgs_ndim, new_bfgs => bfgs, lin_bfgs
USE kinds, ONLY : DP
USE brilz, ONLY : alat, at, bg
@ -43,8 +46,8 @@ SUBROUTINE move_ions()
USE symme, ONLY : s, ftau, nsym, irt
USE ener, ONLY : etot
USE force_mod, ONLY : force
USE control_flags, ONLY : tr2, upscale, iswitch, lbfgs, loldbfgs, &
conv_ions
USE control_flags, ONLY : upscale, iswitch, lbfgs, loldbfgs, lconstrain, &
lmd, conv_ions, alpha0, beta0, tr2
USE relax, ONLY : epse, epsf, starting_scf_threshold
USE cellmd, ONLY : lmovecell, calc
#if defined (__PARA)
@ -52,44 +55,57 @@ SUBROUTINE move_ions()
USE io_global, ONLY : ionode_id
USE mp, ONLY : mp_bcast
#endif
!
! ... external procedures
!
USE constrains_module, ONLY : dist_constrain, check_constrain, &
new_force, compute_penalty
USE basic_algebra_routines, ONLY : norm
!
IMPLICIT NONE
!
REAL(KIND=DP) :: dummy, gv
REAL(KIND=DP), ALLOCATABLE :: dgv(:,:)
REAL(KIND=DP) :: dgv2, theta0
! auxiliar variable :
! gv=0 defines the constrain
! the gradient of gv
! its square modulus
! the value of the one-dimensional const
! ... local variables
!
REAL(KIND=DP), ALLOCATABLE :: tauold(:,:,:)
! previous positions of atoms
REAL(KIND=DP), SAVE :: lambda = 0.5D0
INTEGER :: na
REAL(KIND=DP) :: energy_error, gradient_error
LOGICAL :: step_accepted
LOGICAL :: step_accepted, exst
REAL(KIND=DP), ALLOCATABLE :: pos(:), gradient(:)
!
!
conv_ions = .FALSE.
!
IF ( iswitch == 2 .OR. iswitch == 4 ) THEN
ALLOCATE( tauold( 3, nat, 3 ) )
!
! ... constrains are imposed here
!
IF ( lconstrain ) &
CALL impose_constrains()
!
! ... the file containing old positions is opened ( needed for extrapolation )
!
CALL seqopn( 4, TRIM( prefix ) // '.update', 'FORMATTED', exst )
!
IF ( exst ) THEN
!
ALLOCATE( dgv(3,nat) )
READ( UNIT = 4, FMT = * ) tauold
!
! ... gv is the function which define the constrain, now first of all we
! ... find the constrain theta0 such that gv=0, and find the gradient of
! ... gv, dgv
ELSE
!
dummy = 0.D0
!
CALL constrain( theta0, gv, dgv, dgv2, dummy, nat, tau, alat )
!
! ... find the constrained forces
!
CALL new_force( dgv, dgv2 )
!
DEALLOCATE( dgv )
tauold = 0.D0
!
END IF
!
CLOSE( UNIT = 4, STATUS = 'KEEP' )
!
! ... save the previous two steps ( a total of three steps is saved )
!
tauold(:,:,3) = tauold(:,:,2)
tauold(:,:,2) = tauold(:,:,1)
tauold(:,:,1) = tau(:,:)
!
! ... do the minimization / dynamics step
!
IF ( lmovecell .AND. ( iswitch == 2 .OR. iswitch == 4 ) ) &
@ -107,6 +123,7 @@ SUBROUTINE move_ions()
! ... only one node does the calculation in the parallel case
!
IF ( me == 1 .AND. mypool == 1 ) THEN
!
#endif
!
ALLOCATE( pos( 3 * nat ) )
@ -122,7 +139,7 @@ SUBROUTINE move_ions()
conv_ions )
!
ELSE
!
!
CALL lin_bfgs( pos, etot, gradient, tmp_dir, stdout, epse, &
epsf, energy_error, gradient_error, step_accepted, &
conv_ions )
@ -149,19 +166,20 @@ SUBROUTINE move_ions()
tau = RESHAPE( SOURCE = pos, SHAPE = (/ 3 , nat /) ) / alat
force = - RESHAPE( SOURCE = gradient, SHAPE = (/ 3 , nat /) )
!
CALL output_tau(conv_ions)
CALL output_tau( conv_ions )
!
DEALLOCATE( pos )
DEALLOCATE( gradient )
!
#if defined (__PARA)
!
END IF
!
! ... broadcast calculated quantities to all nodes
!
CALL mp_bcast( tau, ionode_id )
CALL mp_bcast( force, ionode_id )
CALL mp_bcast( tr2, ionode_id )
CALL mp_bcast( tau, ionode_id )
CALL mp_bcast( force, ionode_id )
CALL mp_bcast( tr2, ionode_id )
CALL mp_bcast( conv_ions, ionode_id )
#endif
!
@ -175,191 +193,255 @@ SUBROUTINE move_ions()
!
! ... molecular dynamics schemes are used
!
IF ( iswitch == 3 .OR.iswitch == 4 ) THEN
IF ( lmd ) THEN
!
IF ( calc == ' ' ) CALL dynamics() ! verlet dynamics
IF ( calc /= ' ' ) CALL vcsmd() ! variable cell shape md
!
END IF
!
IF ( iswitch > 4 .OR. iswitch <= 0 ) THEN
!
CALL errore( 'move_ions', 'iswitch value not implemented or wrong', 1 )
!
END IF
! ... check if the new positions satisfy the constrain equation
!
! ... check if the new positions satisfy the constrain equation, in
! ... the CP case this is done inside the routine "cp"
!
IF ( iswitch == 2 .OR. iswitch == 4 ) &
CALL check_constrain( theta0 )
IF ( lconstrain ) CALL check_constrain()
!
! ... before leaving check that the new positions still transform
! ... according to the symmetry of the system.
!
CALL checkallsym( nsym, s, nat, tau, ityp, at, bg, nr1, nr2, nr3, irt, ftau )
!
! ... find the best coefficients for the extrapolation of the potential
!
CALL find_alpha_and_beta( nat, tau, tauold, alpha0, beta0 )
!
#if defined (__PARA)
!
IF ( me == 1 ) CALL poolbcast( 1, alpha0 )
IF ( me == 1 ) CALL poolbcast( 1, beta0 )
!
CALL mp_bcast( alpha0, ionode_id )
CALL mp_bcast( beta0, ionode_id )
!
#endif
!
CALL seqopn( 4, TRIM( prefix ) // '.update', 'FORMATTED', exst )
!
WRITE( UNIT = 4, FMT = * ) tauold
!
CLOSE( UNIT = 4, STATUS = 'KEEP' )
!
DEALLOCATE( tauold )
!
RETURN
!
CONTAINS
!
! ... internal procedures
!
!-----------------------------------------------------------------------
SUBROUTINE impose_constrains()
!-----------------------------------------------------------------------
!
USE constrains_module, ONLY : nconstr
!
IMPLICIT NONE
!
! ... local variables
!
INTEGER :: index, na
REAL(KIND=DP) :: gv
REAL(KIND=DP) :: dgv(3,nat)
REAL(KIND=DP) :: dgv2
! gv = 0 defines the constrain
! the gradient of gv
! its square modulus
!
!
IF ( lbfgs ) THEN
!
! ... BFGS case: a penalty function is used
!
CALL compute_penalty( gv, dgv, dgv2 )
!
etot = etot + lambda * gv**2
!
force(:,:) = force(:,:) - 2.D0 * lambda * gv * dgv(:,:)
!
ELSE IF ( lmd ) THEN
!
! ... molecular dynamics case: lagrange multipliers are used
!
! ... find the constrained forces
!
DO index = 1, nconstr
!
CALL dist_constrain( index, gv, dgv, dgv2 )
!
CALL new_force( dgv, dgv2 )
!
END DO
!
WRITE( stdout, '(/5x,"Constrained forces")')
!
DO na = 1, nat
!
WRITE( stdout, '(3F14.8)') force(:,na)
!
END DO
!
END IF
!
END SUBROUTINE impose_constrains
!
!
!-----------------------------------------------------------------------
SUBROUTINE compute_lambda()
!-----------------------------------------------------------------------
!
USE constrains_module, ONLY : constr_tol
!
IMPLICIT NONE
!
! ... local variables
!
LOGICAL :: ltest
REAL(KIND=DP) :: gv
REAL(KIND=DP) :: dgv(3,nat)
REAL(KIND=DP) :: dgv2
! gv = 0 defines the constrain
! the gradient of gv
! its square modulus
!
!
CALL compute_penalty( gv, dgv, dgv2 )
!
IF ( step_accepted ) THEN
!
lambda_loop: DO
!
IF ( ABS( gv ) > constr_tol ) lambda = lambda * 1.1D0
!
ltest = .TRUE.
!
DO na = 1, nat
!
IF ( 2.D0 * lambda * gv * norm( dgv(:,na) ) > 0.05D0 ) &
ltest = .FALSE.
!
END DO
!
IF ( ltest ) EXIT lambda_loop
!
lambda = lambda * 0.5D0
!
END DO lambda_loop
!
END IF
!
WRITE( stdout, '("LAMBDA = ",F14.10)' ) lambda
WRITE( stdout, '("GV = ",F14.10)' ) gv
WRITE( stdout, '("PENALTY = ",F14.10)' ) lambda * gv**2
!
RETURN
!
END SUBROUTINE compute_lambda
!
!
!-----------------------------------------------------------------------
SUBROUTINE find_alpha_and_beta( nat, tau, tauold, alpha0, beta0 )
!-----------------------------------------------------------------------
!
! ... This routine finds the best coefficients alpha0 and beta0 so that
!
! ... | tau(t+dt) - tau' | is minimum, where
!
! ... tau' = alpha0 * ( tau(t) - tau(t-dt) ) +
! ... beta0 * ( tau(t-dt) - tau(t-2*dt) )
!
USE constants, ONLY : eps8
!
IMPLICIT NONE
!
INTEGER :: nat, na, ipol
REAL(KIND=DP) :: chi, alpha0, beta0, tau(3,nat), tauold(3,nat,3)
REAL(KIND=DP) :: a11, a12, a21, a22, b1, b2, c, det
!
! ... solution of the linear system
!
a11 = 0.D0
a12 = 0.D0
a21 = 0.D0
a22 = 0.D0 + eps8
b1 = 0.D0
b2 = 0.D0
c = 0.D0
!
DO na = 1, nat
!
DO ipol = 1, 3
!
a11 = a11 + ( tauold(ipol,na,1) - tauold(ipol,na,2) )**2
!
a12 = a12 + ( tauold(ipol,na,1) - tauold(ipol,na,2) ) * &
( tauold(ipol,na,2) - tauold(ipol,na,3) )
!
a22 = a22 + ( tauold(ipol,na,2) - tauold(ipol,na,3) )**2
!
b1 = b1 - ( tauold(ipol,na,1) - tau(ipol,na) ) * &
( tauold(ipol,na,1) - tauold(ipol,na,2) )
!
b2 = b2 - ( tauold(ipol,na,1) - tau(ipol,na) ) * &
( tauold(ipol,na,2) - tauold(ipol,na,3) )
!
c = c + ( tauold(ipol,na,1) - tau(ipol,na) )**2
!
END DO
!
END DO
!
a21 = a12
!
det = a11 * a22 - a12 * a21
!
IF ( det < 0.D0 ) CALL errore( 'find_alpha_and_beta', ' det < 0', 1 )
!
! ... case det > 0: a well defined minimum exists
!
IF ( det > 0.D0 ) THEN
!
alpha0 = ( b1 * a22 - b2 * a12 ) / det
beta0 = ( a11 * b2 - a21 * b1 ) / det
!
ELSE
!
! ... case det = 0 : the two increments are linearly dependent,
! ... chose solution with beta = 0
! ... ( discard oldest configuration )
!
alpha0 = 1.D0
beta0 = 0.D0
!
IF ( a11 > 0.D0 ) alpha0 = b1 / a11
!
END IF
!
chi = 0.D0
!
DO na = 1, nat
!
DO ipol = 1, 3
!
chi = chi + ( ( 1.D0 + alpha0 ) * tauold(ipol,na,1) + &
( beta0 - alpha0 ) * tauold(ipol,na,2) - &
beta0 * tauold(ipol, na, 3) - tau(ipol,na) )**2
!
END DO
!
END DO
!
!WRITE( stdout, * ) chi, alpha0, beta0
!
RETURN
!
END SUBROUTINE find_alpha_and_beta
!
END SUBROUTINE move_ions
!
!
!----------------------------------------------------------------------------
SUBROUTINE new_force( dg, dg2 )
!----------------------------------------------------------------------------
!
! find the lagrange multiplier lambda for the problem with one const
!
! force*dg
! lambda = - --------,
! |dg|^2
!
! and redefine the forces:
!
! force = force + lambda*dg
!
! where dg is the gradient of the constraint function
!
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
USE basis, ONLY : nat
USE brilz, ONLY : at, bg
USE force_mod, ONLY : force
USE symme, ONLY : s, nsym, irt
!
IMPLICIT NONE
!
INTEGER :: na, i, ipol
REAL(KIND=DP) :: dg(3, nat), lambda, dg2, sum
REAL(KIND=DP), EXTERNAL :: DDOT
!
!
lambda = 0.D0
!
IF ( dg2 /= 0.D0 ) THEN
!
lambda = - DDOT( 3 * nat, force, 1, dg, 1 ) / dg2
!
CALL DAXPY( 3 * nat, lambda, dg, 1, force, 1)
!
IF ( DDOT( 3 * nat, force, 1, dg, 1 )**2 > 1.D-30 ) THEN
!
CALL errore( 'new_force', 'force is not orthogonal to constrain', - 1 )
PRINT *, DDOT( 3 * nat, force, 1, dg, 1 )**2
!
END IF
!
DO ipol = 1, 3
sum = 0.D0
DO na = 1, nat
sum = sum + force(ipol,na)
END DO
!
! ... impose total force = 0
!
DO na = 1, nat
force(ipol,na) = force(ipol,na) - sum / nat
END DO
END DO
!
! ... resymmetrize (should not be needed, but...)
!
IF ( nsym > 1 ) THEN
!
DO na = 1, nat
CALL trnvect( force(1,na), at, bg, - 1 )
END DO
!
CALL symvect( nat, force, nsym, s, irt )
!
DO na = 1, nat
CALL trnvect( force(1,na), at, bg, 1 )
END DO
!
END IF
!
WRITE( stdout, '(/5x,"Constrained forces")')
!
DO na = 1, nat
WRITE( stdout, '(3F14.8)') ( force(i,na) , i = 1, 3 )
END DO
!
END IF
!
RETURN
!
END SUBROUTINE new_force
!
!
!---------------------------------------------------------------------
SUBROUTINE check_constrain( theta0 )
!---------------------------------------------------------------------
!
! update tau so that the constraint equation g=0 is satisfied,
! use the recursion formula:
!
! g(tau)
! tau' = tau - ------------ * dg(tau)
! |dg(tau)|^2
!
! in normal cases the constraint equation should be always satisfied
! the very first iteration.
!
USE io_global, ONLY : stdout
USE kinds, ONLY : DP
USE constants, ONLY : eps16
USE basis, ONLY : nat, ityp, tau, atm
USE brilz, ONLY : alat
!
IMPLICIT NONE
!
REAL(KIND=DP), ALLOCATABLE :: dg(:,:)
REAL(KIND=DP) :: dg2, g, theta0, dummy
INTEGER :: na, i, j
INTEGER, PARAMETER :: maxiter = 250
!
!
ALLOCATE( dg(3,nat) )
!
CALL constrain( dummy, g, dg, dg2, theta0, nat, tau, alat )
!
WRITE( stdout, '(5X,"G = ",1PE9.2," iteration # ",I3)' ) g, 0
!
DO i = 1, maxiter
!
! ... check if g=0
!
IF ( ABS( g ) < eps16 ) GO TO 14
!
! ... if g<>0 find new tau = tau - g*dg/dg2 and check again
!
CALL DAXPY( 3 * nat, - g / dg2, dg, 1, tau, 1 )
!
CALL constrain( dummy, g, dg, dg2, theta0, nat, tau, alat )
!
WRITE( stdout, '(5X,"G = ",1PE9.2," iteration # ",I3)' ) g, i
!
END DO
!
CALL errore( 'new_dtau', 'g=0 is not satisfied g=', - 1 )
!
14 CONTINUE
!
WRITE( stdout, '(5X,"Number of step(s): ",I3)') i - 1
!
! ... if the atomic positions have been corrected write them on output
!
IF ( i > 1 ) THEN
!
WRITE( stdout, '(/5X,"Corrected atomic positions:",/)')
DO na = 1, nat
WRITE( stdout,'(A3,3X,3F14.9)') atm(ityp(na)), ( tau(j,na), j = 1, 3 )
END DO
!
END IF
!
DEALLOCATE( dg )
!
RETURN
!
END SUBROUTINE check_constrain

19
PW/update_pot.f90
View File

@ -45,12 +45,18 @@ SUBROUTINE update_pot()
!
CALL start_clock( 'update_pot' )
!
IF ( order == 0 ) RETURN
IF ( order == 0 ) THEN
!
CALL stop_clock( 'update_pot' )
!
RETURN
!
END IF
!
IF ( order > 2 .AND. ( lbfgs .OR. lneb ) ) THEN
IF ( order > 2 .AND. lneb ) THEN
!
order = 2
CALL errore( 'update_pot', 'order > 2 not allowed in bfgs', -1 )
CALL errore( 'update_pot', 'order > 2 not allowed in neb', -1 )
!
END IF
!
@ -327,7 +333,7 @@ SUBROUTINE extrapolate_wfcs()
!
! ... extrapolate the wfc's,
! ... if order=3 use the second order extrapolation formula
! ... alpha0 and beta0 are calculated in "dynamics"
! ... alpha0 and beta0 are calculated in "move_ions"
!
IF ( order > 2 ) THEN
!
@ -342,19 +348,22 @@ SUBROUTINE extrapolate_wfcs()
CALL davcio( evcold, nwordwfc, iunoldwfc2, ik, - 1 )
!
evc = evc - beta0 * evcold
!
END IF
!
ELSE
!
evc = 2 * evcold - evc
evc = 2.D0 * evcold - evc
!
END IF
!
! ... move the files: "old" -> "old1" and "now" -> "old"
!
IF ( order > 2 ) THEN
!
CALL davcio( evcold, nwordwfc, iunoldwfc, ik, - 1 )
CALL davcio( evcold, nwordwfc, iunoldwfc2, ik, 1 )
!
END IF
!
CALL davcio( evcold, nwordwfc, iunwfc, ik, - 1 )

871
PW/vcsmd.f90
View File

@ -1,14 +1,16 @@
!
! Copyright (C) 2001 PWSCF group
! Copyright (C) 2001-2004 PWSCF group
! 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 .
!
!-----------------------------------------------------------------------
subroutine vcsmd
!-----------------------------------------------------------------------
#include "machine.h"
!
!----------------------------------------------------------------------------
SUBROUTINE vcsmd()
!----------------------------------------------------------------------------
!
! Main (interface) routine between PWSCF and the variable-cell shape
! molecular dynamics code by R.M. Wentzcovitch, PRB 44, 2358 (1991).
!
@ -23,26 +25,30 @@ subroutine vcsmd
! calc = 'nm' : Wentzcovitch's new cell minimization by damped dynam
!
! Dynamics performed using Beeman algorithm, J. Comp. Phys. 20, 130 (1976))
!
#include "machine.h"
USE io_global, ONLY : stdout
use constants, only : e2, uakbar
use brilz
use basis
use cellmd
use dynam
use relax
use force_mod
USE control_flags
use ener, only: etot
use io_files, only : prefix
#ifdef __PARA
use para
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE constants, ONLY : e2, uakbar
USE brilz, ONLY : omega, alat, at, bg
USE basis, ONLY : tau, nat, ntyp, ityp, atm
USE cellmd, ONLY : nzero, ntimes, calc, press, at_old, omega_old, &
cmass, ttol, ntcheck, lmovecell
USE dynam, ONLY : dt, temperature, amass
USE relax, ONLY : epse, epsf, fixatom
USE force_mod, ONLY : force, sigma
USE control_flags, ONLY : istep, conv_ions
USE parameters, ONLY : ntypx
USE ener, ONLY : etot
USE io_files, ONLY : prefix
USE parser, ONLY : delete_if_present
#if defined (__PARA)
USE para, ONLY : me, mypool
#endif
implicit none
!
! I/O variable first
IMPLICIT NONE
!
! ... I/O variable first
!
! PWSCF variables
! nat = total number of atoms
@ -55,548 +61,383 @@ subroutine vcsmd
! cmass = cell mass in ryd units.
! press = target pressure in ryd/(a.u.)^3
!
! local variable
! ... local variables
!
real(kind=DP) :: p, & ! virial pressure
REAL(KIND=DP) :: p, & ! virial pressure
vcell, & ! cell volume
avec (3, 3), & ! at(3,3) * alat
aveci (3, 3), & ! avec at t-dt
avecd (3, 3), & ! d(avec)/dt
avec2d (3, 3),& ! d2(avec)/dt2
avec2di (3, 3),&! d2(avec)/dt2 at t-dt
avec0 (3, 3), & ! avec at t = 0
sig0 (3, 3), & ! sigma at t=0
avec(3,3), & ! at(3,3) * alat
aveci(3,3), & ! avec at t-dt
avecd(3,3), & ! d(avec)/dt
avec2d(3,3), & ! d2(avec)/dt2
avec2di(3,3), & ! d2(avec)/dt2 at t-dt
avec0(3,3), & ! avec at t = 0
sig0(3,3), & ! sigma at t=0
v0 ! volume at t=0
real(kind=DP), allocatable :: &
rat (:, :), & ! atomic positions (lattice coord)
rati (:, :), & ! rat at previous step
ratd (:, :), & ! rat derivatives at current step
rat2d (:, :), & ! rat 2nd derivatives at current step
rat2di (:,:), & ! rat 2nd derivatives at previous step
tauold (:, :, :)! additional history variables
real(kind=DP) :: &
avmod (3), theta (3, 3), & ! used to monitor cell dynamics
enew, e_start, & ! DFT energy at current and first step
eold, & ! DFT energy at previous step
REAL(KIND=DP), ALLOCATABLE :: &
rat(:,:), & ! atomic positions (lattice coord)
rati(:,:), & ! rat at previous step
ratd(:,:), & ! rat derivatives at current step
rat2d(:,:), & ! rat 2nd derivatives at current step
rat2di(:,:), & ! rat 2nd derivatives at previous step
tauold(:,:,:) ! additional history variables
REAL(KIND=DP) :: &
avmod(3), theta(3,3), & ! used to monitor cell dynamics
enew, e_start, & ! DFT energy at current and first step
eold, & ! DFT energy at previous step
uta, eka, eta, ekla, utl, etl, ut, ekint, edyn, & ! other energies
acu, ack, acp, acpv, avu, avk, avp, avpv, & ! acc.& avrg. ener
tnew, pv, & ! istantaneous temperature and p*vcell
sigmamet (3, 3), & ! sigma = avec^-1 * vcell = bg/alat*omega
vx2 (ntypx), vy2 (ntypx), vz2 (ntypx), & ! work vectors
vmean (ntypx), rms (ntypx), ekin (ntypx), & ! work vectors
tnew, pv, & ! istantaneous temperature and p*vcell
sigmamet(3,3), & ! sigma = avec^-1 * vcell = bg/alat*omega
vx2(ntypx), vy2(ntypx), vz2(ntypx), & ! work vectors
vmean(ntypx), rms(ntypx), ekin(ntypx), & ! work vectors
tempo, time_au, epsp
character(len=3) :: ios ! status (old or new) for I/O files
character(len=6) :: ipos ! status ('append' or 'asis') for I/O files
logical :: exst
integer :: idone, na, nst, ipol, i, j, k
! last completed time step
! counter on atoms
! counter on completed moves
CHARACTER(LEN=3) :: ios ! status (old or new) for I/O files
CHARACTER(LEN=6) :: ipos ! status ('append' or 'asis') for I/O files
LOGICAL :: exst
INTEGER :: idone, na, nst, ipol, i, j, k
! last completed time step
! counter on atoms
! counter on completed moves
!
! I/O units
! ... I/O units
!
integer :: iun_e, iun_eal, iun_ave, iun_p, iun_avec, iun_tv
parameter (iun_e=21, iun_eal=22, iun_ave=23, iun_p=24, iun_avec=25, iun_tv=26)
INTEGER, PARAMETER :: iun_e = 21, &
iun_eal = 22, &
iun_ave = 23, &
iun_p = 24, &
iun_avec = 25, &
iun_tv = 26
!
! Allocate work arrays
!
allocate (rat(3,nat), rati(3,nat), ratd(3,nat), rat2d(3,nat), rat2di(3,nat),&
tauold(3,nat,3))
! ... Allocate work arrays
!
! open MD history file (if not present this is a new run!)
ALLOCATE( rat(3,nat) )
ALLOCATE( rati(3,nat) )
ALLOCATE( ratd(3,nat) )
ALLOCATE( rat2d(3,nat) )
ALLOCATE( rat2di(3,nat) )
ALLOCATE( tauold(3,nat,3) )
!
call seqopn (4, trim(prefix)//'.md', 'formatted', exst)
if (.not.exst) then
close (unit = 4, status = 'delete')
if (istep.ne.1) call errore ('vcsmd', ' previous MD history got lost', 1)
tnew = 0.d0
acu = 0.d0
ack = 0.d0
acp = 0.d0
acpv = 0.d0
avu = 0.d0
avk = 0.d0
avp = 0.d0
avpv = 0.d0
! ... open MD history file (if not present this is a new run!)
!
CALL seqopn( 4, TRIM( prefix ) // '.md', 'FORMATTED', exst )
!
IF ( .NOT. exst ) THEN
!
CLOSE( UNIT = 4, STATUS = 'DELETE' )
IF ( istep /= 1 ) &
CALL errore( 'vcsmd', ' previous MD history got lost', 1 )
!
tnew = 0.D0
acu = 0.D0
ack = 0.D0
acp = 0.D0
acpv = 0.D0
avu = 0.D0
avk = 0.D0
avp = 0.D0
avpv = 0.D0
nzero = 0
tauold(:,:,:) = 0.d0
eold = etot + 2*epse ! set value for eold at first iteration
else
tauold(:,:,:) = 0.D0
!
! read MD run history (idone is the last completed MD step)
! ... set value for eold at first iteration
!
read (4, * ) rati, ratd, rat2d, rat2di, tauold
read (4, * ) aveci, avecd, avec2d, avec2di
read (4, * ) avec0, sig0, v0, e_start, eold
read (4, * ) acu, ack, acp, acpv, avu, avk, avp, avpv, sigmamet
read (4, * ) idone, nzero, ntimes
close (unit = 4, status = 'keep')
istep = idone+1
call DCOPY (3 * nat, tauold (1, 1, 2), 1, tauold (1, 1, 3), 1)
call DCOPY (3 * nat, tauold (1, 1, 1), 1, tauold (1, 1, 2), 1)
endif
eold = etot + 2.D0 * epse
!
ELSE
!
! ... read MD run history (idone is the last completed MD step)
!
READ( 4, * ) rati, ratd, rat2d, rat2di, tauold
READ( 4, * ) aveci, avecd, avec2d, avec2di
READ( 4, * ) avec0, sig0, v0, e_start, eold
READ( 4, * ) acu, ack, acp, acpv, avu, avk, avp, avpv, sigmamet
READ( 4, * ) idone, nzero, ntimes
!
CLOSE( UNIT = 4, STATUS = 'KEEP' )
!
istep = idone + 1
!
tauold(:,:,3) = tauold(:,:,2)
tauold(:,:,2) = tauold(:,:,1)
!
END IF
!
! check if convergence for structural minimization is achieved
! ... check if convergence for structural minimization is achieved
!
if (calc.eq.'mm') then
conv_ions = eold - etot .lt. epse
do i = 1, 3
do na = 1, nat - fixatom
conv_ions = conv_ions.and.abs (force (i, na) ) .lt. epsf
enddo
enddo
if (conv_ions) then
WRITE( stdout,'(/5x,"Damped Dynamics: convergence achieved, Efinal=",&
& f15.8)') etot
call output_tau(.TRUE.)
return
end if
end if
if (calc.eq.'nm' .or. calc.eq.'cm') then
epsp = 0.5 ! kbar
conv_ions = eold - etot .lt. epse
do i = 1, 3
do na = 1, nat - fixatom
conv_ions = conv_ions.and.abs (force (i, na) ) .lt. epsf
enddo
enddo
do i =1,3
conv_ions = conv_ions .and. abs (sigma(i,i) - press)*uakbar.lt. epsp
do j =i+1,3
conv_ions = conv_ions .and. abs ( sigma(i,j) )*uakbar.lt. epsp
end do
end do
if (conv_ions) then
if (calc.eq.'cm') WRITE( stdout, &
'(/5x,"Parrinello-Rahman Damped Dynamics: convergence achieved, ", &
& "Efinal=", f15.8)') etot
if (calc.eq.'nm') WRITE( stdout, &
'(/5x,"Wentzcovitch Damped Dynamics: convergence achieved, ", &
& "Efinal=", f15.8)') etot
WRITE( stdout,'(/72("-")//5x,"Final estimate of lattice vectors ", &
& "(input alat units)")')
WRITE( stdout, '(3f14.9)') ( (at (i, k) , i = 1, 3) , k = 1, 3)
WRITE( stdout,'(" final unit-cell volume =",f12.4," (a.u.)^3")') omega
WRITE( stdout,'(" input alat = ",f12.4," (a.u.)")') alat
IF ( calc == 'mm' ) THEN
!
conv_ions = ( eold - etot ) < epse
!
DO i = 1, 3
DO na = 1, nat
conv_ions = conv_ions .AND. ( ABS( force(i,na) ) < epsf )
END DO
END DO
!
IF ( conv_ions ) THEN
!
call output_tau (.TRUE.)
WRITE( UNIT = stdout, &
FMT = '(/,5X,"Damped Dynamics: convergence achieved, Efinal=", &
&F15.8)' ) etot
!
return
end if
end if
call DCOPY (3 * nat, tau, 1, tauold (1, 1, 1), 1)
CALL output_tau( .TRUE. )
!
RETURN
!
END IF
!
ELSE IF ( calc == 'nm' .OR. calc == 'cm' ) THEN
!
epsp = 0.5D0 ! kbar
!
conv_ions = ( eold - etot ) < epse
!
DO i = 1, 3
DO na = 1, nat
conv_ions = conv_ions .AND. ( ABS( force(i,na) ) < epsf )
END DO
END DO
!
DO i = 1, 3
!
conv_ions = conv_ions .AND. &
( ABS( sigma(i,i) - press) * uakbar < epsp )
!
DO j = ( i + 1 ), 3
conv_ions = conv_ions .AND. &
( ABS( sigma(i,j) ) * uakbar < epsp )
END DO
!
END DO
!
IF ( conv_ions ) THEN
!
IF ( calc == 'cm' ) &
WRITE( UNIT = stdout, &
FMT = '(/5X,"Parrinello-Rahman Damped Dynamics: ", &
& "convergence achieved, Efinal=", F15.8)' ) etot
IF ( calc == 'nm' ) &
WRITE( UNIT = stdout, &
FMT = '(/5X,"Wentzcovitch Damped Dynamics: ", &
& "convergence achieved, Efinal=", F15.8)' ) etot
!
WRITE( UNIT = stdout, &
FMT = '(/72("-")//5X,"Final estimate of lattice vectors ", &
& "(input alat units)")' )
WRITE( UNIT = stdout, &
FMT = '(3F14.9)') ( ( at(i,k) , i = 1, 3 ) , k = 1, 3 )
WRITE( UNIT = stdout, &
FMT = '(" final unit-cell volume =",F12.4," (a.u.)^3")') omega
WRITE( UNIT = stdout, &
FMT = '(" input alat = ",F12.4," (a.u.)")') alat
!
CALL output_tau( .TRUE. )
!
RETURN
!
END IF
!
END IF
!
tauold(:,:,1) = tau(:,:)
!
time_au = 0.0000242 * e2
tempo = (istep - 1) * dt * time_au
if (istep.eq.1 .and. calc.eq.'mm') &
WRITE( stdout,'(/5x,"Damped Dynamics Minimization", /5x, &
& "convergence thresholds: EPSE = ", e8.2," EPSF = ",e8.2)') &
!
tempo = ( istep - 1 ) * dt * time_au
!
IF ( istep == 1 .AND. calc == 'mm' ) &
WRITE( stdout,'(/5X,"Damped Dynamics Minimization",/5X, &
& "convergence thresholds: EPSE = ",E8.2," EPSF = ",E8.2)' ) &
epse, epsf
if (istep.eq.1 .and. calc.eq.'cm') &
WRITE( stdout,'(/5x,"Parrinello-Rahman Damped Cell-Dynamics Minimization", &
& /5x, "convergence thresholds: EPSE = ", e8.2," EPSF = ", &
& e8.2," EPSP = ",e8.2 )') epse, epsf, epsp
if (istep.eq.1 .and. calc.eq.'nm') &
WRITE( stdout,'(/5x,"Wentzcovitch Damped Cell-Dynamics Minimization", /5x, &
& "convergence thresholds: EPSE = ", e8.2," EPSF = ",e8.2,&
& " EPSP = ",e8.2 )') epse, epsf, epsp
WRITE( stdout, '(/5x,"Entering Dynamics; it = ",i5," time = ", &
& f8.5," pico-seconds"/)') istep, tempo
! WRITE( stdout,*) ' enter vcsmd ', istep,idone,exst
IF ( istep == 1 .AND. calc == 'cm' ) &
WRITE( stdout,'(/5X,"Parrinello-Rahman Damped Cell-Dynamics Minimization", &
& /5X, "convergence thresholds: EPSE = ",E8.2," EPSF = ", &
& E8.2," EPSP = ",E8.2 )' ) epse, epsf, epsp
IF ( istep == 1 .AND. calc == 'nm' ) &
WRITE( stdout,'(/5X,"Wentzcovitch Damped Cell-Dynamics Minimization", /5x, &
& "convergence thresholds: EPSE = ",E8.2," EPSF = ",E8.2, &
& " EPSP = ",E8.2 )' ) epse, epsf, epsp
!
! save cell shape of previous step
WRITE( stdout, '(/5X,"Entering Dynamics; it = ",I5," time = ", &
& F8.5," pico-seconds"/)' ) istep, tempo
!
! ... save cell shape of previous step
!
at_old = at
!
call DCOPY (9, at, 1, at_old, 1)
omega_old = omega
!
! Translate
! ... Translate
!
! define rat as the atomic positions in lattice coordinates
! ... define rat as the atomic positions in lattice coordinates
!
call DCOPY (3 * nat, tau, 1, rat, 1)
call cryst_to_cart (nat, rat, bg, - 1)
rat = tau
!
call DCOPY (9, at, 1, avec, 1)
call DSCAL (9, alat, avec, 1)
CALL cryst_to_cart( nat, rat, bg, -1 )
!
! convert forces to lattice coordinates
avec = alat * at
!
call cryst_to_cart (nat, force, bg, - 1)
call DSCAL (3 * nat, 1.d0 / alat, force, 1)
! ... convert forces to lattice coordinates
!
! scale stress to stress*omega
CALL cryst_to_cart( nat, force, bg, -1 )
!
force = force / alat
!
! ... scale stress to stress*omega
!
sigma = sigma * omega
!
call DSCAL (9, omega, sigma, 1)
vcell = omega
if (istep.eq.1) then
!
IF ( istep == 1 ) THEN
!
e_start = etot
!
enew = etot - e_start
#ifdef DEBUG_VCSMD
write (45,*) 'istep=',istep
write (45,*) 'ntyp=', ntyp
write (45,*) 'nat=', nat
write (45,*) 'rat=',rat
write (45,*) 'ityp=',ityp
write (45,*) 'avec=',avec
write (45,*) 'vcell=',vcell
write (45,*) 'force=',force
write (45,*) 'sigma=',sigma
write (45,*) 'calc=',calc
write (45,*) 'temperature=',temperature
write (45,*) 'vx2=',vx2
write (45,*) 'vy2=',vy2
write (45,*) 'vz2=',vz2
write (45,*) 'rms=',rms
write (45,*) 'vmean=',vmean
write (45,*) 'ekin=',ekin
write (45,*) 'avmod=',avmod
write (45,*) 'theta=',theta
write (45,*) 'amass=',amass
write (45,*) 'cmass=',cmass
write (45,*) 'press=',press
write (45,*) 'p=',p
write (45,*) 'dt=',dt
write (45,*) 'aveci=',aveci
write (45,*) 'avecd=',avecd
write (45,*) 'avec2d=',avec2d
write (45,*) 'avec2di=',avec2di
write (45,*) 'sigmamet=',sigmamet
write (45,*) 'sig0=', sig0
write (45,*) 'avec0=',avec0
write (45,*) 'v0=',v0
write (45,*) 'rati=',rati
write (45,*) 'ratd=',ratd
write (45,*) 'rat2d=',rat2d
write (45,*) 'rat2di=',rat2di
write (45,*) 'enew=',enew
write (45,*) 'uta=',uta
write (45,*) 'eka=',eka
write (45,*) 'eta=',eta
write (45,*) 'ekla=',ekla
write (45,*) 'utl=',utl
write (45,*) 'etl=',etl
write (45,*) 'ut=',ut
write (45,*) 'ekint=',ekint
write (45,*) 'edyn=',edyn
#endif
call init (ntyp, nat, ntyp, nat-fixatom, rat, ityp, avec, vcell, force, &
sigma, calc, temperature, vx2, vy2, vz2, rms, vmean, ekin, &
avmod, theta, amass, cmass, press, p, dt, aveci, avecd, avec2d, &
avec2di, sigmamet, sig0, avec0, v0, rati, ratd, rat2d, rat2di, &
enew, uta, eka, eta, ekla, utl, etl, ut, ekint, edyn)
#ifdef DEBUG_VCSMD
write (46,*) 'istep=',istep
write (46,*) 'ntyp=', ntyp
write (46,*) 'nat=', nat
write (46,*) 'rat=',rat
write (46,*) 'ityp=',ityp
write (46,*) 'avec=',avec
write (46,*) 'vcell=',vcell
write (46,*) 'force=',force
write (46,*) 'sigma=',sigma
write (46,*) 'calc=',calc
write (46,*) 'temperature=',temperature
write (46,*) 'vx2=',vx2
write (46,*) 'vy2=',vy2
write (46,*) 'vz2=',vz2
write (46,*) 'rms=',rms
write (46,*) 'vmean=',vmean
write (46,*) 'ekin=',ekin
write (46,*) 'avmod=',avmod
write (46,*) 'theta=',theta
write (46,*) 'amass=',amass
write (46,*) 'cmass=',cmass
write (46,*) 'press=',press
write (46,*) 'p=',p
write (46,*) 'dt=',dt
write (46,*) 'aveci=',aveci
write (46,*) 'avecd=',avecd
write (46,*) 'avec2d=',avec2d
write (46,*) 'avec2di=',avec2di
write (46,*) 'sigmamet=',sigmamet
write (46,*) 'sig0=', sig0
write (46,*) 'avec0=',avec0
write (46,*) 'v0=',v0
write (46,*) 'rati=',rati
write (46,*) 'ratd=',ratd
write (46,*) 'rat2d=',rat2d
write (46,*) 'rat2di=',rat2di
write (46,*) 'enew=',enew
write (46,*) 'uta=',uta
write (46,*) 'eka=',eka
write (46,*) 'eta=',eta
write (46,*) 'ekla=',ekla
write (46,*) 'utl=',utl
write (46,*) 'etl=',etl
write (46,*) 'ut=',ut
write (46,*) 'ekint=',ekint
write (46,*) 'edyn=',edyn
#endif
else
! WRITE( stdout,'(3f12.6)') ((ratd(ipol,na),ipol=1,3),na=1,nat)
!
CALL init( ntyp, nat, ntyp, nat-fixatom, rat, ityp, avec, vcell, force, &
sigma, calc, temperature, vx2, vy2, vz2, rms, vmean, ekin, &
avmod, theta, amass, cmass, press, p, dt, aveci, avecd, avec2d,&
avec2di, sigmamet, sig0, avec0, v0, rati, ratd, rat2d, rat2di, &
enew, uta, eka, eta, ekla, utl, etl, ut, ekint, edyn )
!
ELSE
!
enew = etot - e_start
#ifdef DEBUG_VCSMD
write (45,*) 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
write (45,*) 'istep=', istep
write (45,*) 'ntyp=',ntyp
write (45,*) 'nat=',nat
write (45,*) 'ityp=',ityp
write (45,*) 'rat=',rat
write (45,*) 'avec=',avec
write (45,*) 'vcell=',vcell
write (45,*) 'force=',force
write (45,*) 'sigma=',sigma
write (45,*) 'calc=',calc
write (45,*) 'avmod=',avmod
write (45,*) 'theta=',theta
write (45,*) 'amass=',amass
write (45,*) 'cmass=',cmass
write (45,*) 'press=',press
write (45,*) 'p=',p
write (45,*) 'dt=',dt
write (45,*) 'avecd=',avecd
write (45,*) 'avec2d=',avec2d
write (45,*) 'aveci=',aveci
write (45,*) 'avec2di=',avec2di
write (45,*) 'sigmamet=',sigmamet
write (45,*) 'sig0=',sig0
write (45,*) 'avec0=',avec0
write (45,*) 'v0=',v0
write (45,*) 'ratd=',ratd
write (45,*) 'rat2d=',rat2d
write (45,*) 'rati=',rati
write (45,*) 'rat2di=',rat2di
write (45,*) 'enew=',enew
write (45,*) 'uta=',uta
write (45,*) 'eka=',eka
write (45,*) 'eta=',eta
write (45,*) 'ekla=',ekla
write (45,*) 'utl=',utl
write (45,*) 'etl=',etl
write (45,*) 'ut=',ut
write (45,*) 'ekint=',ekint
write (45,*) 'edyn=',edyn
write (45,*) 'temperature=',temperature
write (45,*) 'ttol=',ttol
write (45,*) 'ntcheck=',ntcheck
write (45,*) 'ntimes=',ntimes
write (45,*) 'istep=',istep
write (45,*) 'tnew=',tnew
write (45,*) 'nzero=',nzero
write (45,*) 'nat=',nat
write (45,*) 'acu=',acu
write (45,*) 'ack=',ack
write (45,*) 'acp=',acp
write (45,*) 'acpv=',acpv
write (45,*) 'avu=',avu
write (45,*) 'avk=',avk
write (45,*) 'avp=',avp
write (45,*) 'avpv=',avpv
#endif
call move (ntyp, nat, ntyp, ityp, rat, avec, vcell, force, &
sigma, calc, avmod, theta, amass, cmass, press, p, dt, avecd, &
avec2d, aveci, avec2di, sigmamet, sig0, avec0, v0, ratd, rat2d, &
rati, rat2di, enew, uta, eka, eta, ekla, utl, etl, ut, ekint, &
edyn, temperature, ttol, ntcheck, ntimes, istep, tnew, nzero, &
nat-fixatom, acu, ack, acp, acpv, avu, avk, avp, avpv)
#ifdef DEBUG_VCSMD
write (46,*) 'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'
write (46,*) 'istep=', istep
write (46,*) 'ntyp=',ntyp
write (46,*) 'nat=',nat
write (46,*) 'ityp=',ityp
write (46,*) 'rat=',rat
write (46,*) 'avec=',avec
write (46,*) 'vcell=',vcell
write (46,*) 'force=',force
write (46,*) 'sigma=',sigma
write (46,*) 'calc=',calc
write (46,*) 'avmod=',avmod
write (46,*) 'theta=',theta
write (46,*) 'amass=',amass
write (46,*) 'cmass=',cmass
write (46,*) 'press=',press
write (46,*) 'p=',p
write (46,*) 'dt=',dt
write (46,*) 'avecd=',avecd
write (46,*) 'avec2d=',avec2d
write (46,*) 'aveci=',aveci
write (46,*) 'avec2di=',avec2di
write (46,*) 'sigmamet=',sigmamet
write (46,*) 'sig0=',sig0
write (46,*) 'avec0=',avec0
write (46,*) 'v0=',v0
write (46,*) 'ratd=',ratd
write (46,*) 'rat2d=',rat2d
write (46,*) 'rati=',rati
write (46,*) 'rat2di=',rat2di
write (46,*) 'enew=',enew
write (46,*) 'uta=',uta
write (46,*) 'eka=',eka
write (46,*) 'eta=',eta
write (46,*) 'ekla=',ekla
write (46,*) 'utl=',utl
write (46,*) 'etl=',etl
write (46,*) 'ut=',ut
write (46,*) 'ekint=',ekint
write (46,*) 'edyn=',edyn
write (46,*) 'temperature=',temperature
write (46,*) 'ttol=',ttol
write (46,*) 'ntcheck=',ntcheck
write (46,*) 'ntimes=',ntimes
write (46,*) 'istep=',istep
write (46,*) 'tnew=',tnew
write (46,*) 'nzero=',nzero
write (46,*) 'nat=',nat
write (46,*) 'acu=',acu
write (46,*) 'ack=',ack
write (46,*) 'acp=',acp
write (46,*) 'acpv=',acpv
write (46,*) 'avu=',avu
write (46,*) 'avk=',avk
write (46,*) 'avp=',avp
write (46,*) 'avpv=',avpv
#endif
endif
!
CALL move( ntyp, nat, ntyp, ityp, rat, avec, vcell, force, &
sigma, calc, avmod, theta, amass, cmass, press, p, dt, avecd, &
avec2d, aveci, avec2di, sigmamet, sig0, avec0, v0, ratd, rat2d,&
rati, rat2di, enew, uta, eka, eta, ekla, utl, etl, ut, ekint, &
edyn, temperature, ttol, ntcheck, ntimes, istep, tnew, nzero, &
nat-fixatom, acu, ack, acp, acpv, avu, avk, avp, avpv)
!
END IF
!
pv = p * omega
!
! write on output files several control quantities
!
#ifdef __PARA
!
! only the first processor does this I/O
!
if (me.eq.1.and.mypool.eq.1) then
IF ( calc /= 'mm' .AND. calc /= 'nm' .AND. calc /= 'cm' ) THEN
!
! ... write on output files several control quantities
!
#if defined (__PARA)
!
! ... only the first processor does this I/O
!
IF ( me == 1 .AND. mypool == 1 ) THEN
#endif
!
! NB: at the first iteration files should not be present,
! for subsequent iterations they should.
!
if (istep.eq.1) then
call delete_if_present ( 'e')
call delete_if_present ( 'eal')
call delete_if_present ( 'ave')
call delete_if_present ( 'p')
call delete_if_present ( 'avec')
call delete_if_present ( 'tv')
ios = 'new'
ipos= 'asis'
else
ios = 'old'
ipos= 'append'
endif
open (unit=iun_e, file='e', status=ios,form='formatted',position=ipos)
open (unit=iun_eal, file='eal', status=ios,form='formatted',position=ipos)
open (unit=iun_ave, file='ave', status=ios,form='formatted',position=ipos)
open (unit=iun_p, file='p', status=ios,form='formatted',position=ipos)
open (unit=iun_avec,file='avec',status=ios,form='formatted',position=ipos)
open (unit=iun_tv, file='tv', status=ios,form='formatted',position=ipos)
nst = istep - 1
write (iun_e, 101) ut, ekint, edyn, pv, nst
write (iun_eal, 103) uta, eka, eta, utl, ekla, etl, nst
write (iun_ave, 104) avu, avk, nst
write (iun_p, 105) press, p, avp, nst
if (calc (1:1) .ne.'m') write (iun_avec, 103) (avmod(k), k=1,3), &
theta(1,2), theta(2,3), theta(3,1), nst
write (iun_tv, 104) vcell, tnew, nst
close (unit = iun_e, status = 'keep')
close (unit = iun_eal, status = 'keep')
close (unit = iun_ave, status = 'keep')
close (unit = iun_p, status = 'keep')
close (unit = iun_avec, status = 'keep')
close (unit = iun_tv, status = 'keep')
#ifdef __PARA
endif
!
! ... NB: at the first iteration files should not be present,
! ... for subsequent iterations they should.
!
IF ( istep == 1 ) THEN
!
CALL delete_if_present( 'e' )
CALL delete_if_present( 'eal' )
CALL delete_if_present( 'ave' )
CALL delete_if_present( 'p' )
CALL delete_if_present( 'avec' )
CALL delete_if_present( 'tv' )
!
ios = 'NEW'
ipos = 'ASIS'
!
ELSE
!
ios = 'OLD'
ipos = 'APPEND'
!
END IF
!
OPEN( UNIT = iun_e, FILE = 'e', STATUS = ios, &
FORM = 'FORMATTED', POSITION = ipos )
OPEN( UNIT = iun_eal, FILE = 'eal', STATUS = ios, &
FORM = 'FORMATTED', POSITION = ipos )
OPEN( UNIT = iun_ave, FILE = 'ave', STATUS = ios, &
FORM = 'FORMATTED', POSITION = ipos )
OPEN( UNIT = iun_p, FILE = 'p', STATUS = ios, &
FORM = 'FORMATTED', POSITION = ipos )
OPEN( UNIT = iun_avec, FILE = 'avec', STATUS = ios, &
FORM = 'FORMATTED', POSITION = ipos )
OPEN( UNIT = iun_tv, FILE = 'tv', STATUS = ios, &
FORM = 'FORMATTED', POSITION = ipos )
!
nst = istep - 1
!
WRITE( iun_e, 101 ) ut, ekint, edyn, pv, nst
WRITE( iun_eal, 103 ) uta, eka, eta, utl, ekla, etl, nst
WRITE( iun_ave, 104 ) avu, avk, nst
WRITE( iun_p, 105 ) press, p, avp, nst
!
IF ( calc(1:1) /= 'm' ) &
WRITE( iun_avec, 103 ) &
avmod(:), theta(1,2), theta(2,3), theta(3,1), nst
!
WRITE( iun_tv, 104 ) vcell, tnew, nst
!
CLOSE( UNIT = iun_e, STATUS = 'KEEP' )
CLOSE( UNIT = iun_eal, STATUS = 'KEEP' )
CLOSE( UNIT = iun_ave, STATUS = 'KEEP' )
CLOSE( UNIT = iun_p, STATUS = 'KEEP' )
CLOSE( UNIT = iun_avec, STATUS = 'KEEP' )
CLOSE( UNIT = iun_tv, STATUS = 'KEEP' )
!
#if defined (__PARA)
END IF
#endif
!
! update configuration in PWSCF variables
END IF
!
call DCOPY (9, avec, 1, at, 1)
call DSCAL (9, 1.d0 / alat, at, 1)
call volume (alat, at (1, 1), at (1, 2), at (1, 3), omega)
call recips (at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2) , bg(1,3) )
call DCOPY (3 * nat, rat, 1, tau, 1)
if (lmovecell) then
! ... update configuration in PWSCF variables
!
at = avec / alat
!
CALL volume( alat, at(1,1), at(1,2), at(1,3), omega )
!
CALL recips( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2) , bg(1,3) )
!
tau = rat
!
IF ( lmovecell ) THEN
!
WRITE( stdout, * ) ' new lattice vectors (alat unit) :'
WRITE( stdout, '(3f14.9)') ( (at (i, k) , i = 1, 3) , k = 1, 3)
WRITE( stdout,'(a,f12.4,a)') ' new unit-cell volume =', omega, ' (a.u.)^3'
endif
WRITE( stdout, '(3F14.9)') ( ( at(i,k) , i = 1, 3 ) , k = 1, 3 )
WRITE( stdout,'(A,F12.4,A)') ' new unit-cell volume =', omega, ' (a.u.)^3'
!
END IF
!
WRITE( stdout, * ) ' new positions in cryst coord'
WRITE( stdout,'(a3,3x,3f14.9)') (atm(ityp(na)), (tau(i,na), i=1,3),na=1,nat)
WRITE( stdout,'(A3,3X,3F14.9)') ( atm(ityp(na)), tau(:,na), na = 1, nat )
WRITE( stdout, * ) ' new positions in cart coord (alat unit)'
call cryst_to_cart (nat, tau, at, + 1)
WRITE( stdout,'(a3,3x,3f14.9)') (atm(ityp(na)), (tau(i,na), i=1,3),na=1,nat)
WRITE( stdout, '(/5x,"Ekin = ",f14.8," Ryd T = ",f6.1," K ", &
& " Etot = ",f14.8)') ekint, tnew, edyn + e_start
!
! save MD history on file
CALL cryst_to_cart( nat, tau, at, 1 )
!
call seqopn (4, trim(prefix)//'.md', 'formatted', exst)
write (4, * ) rati, ratd, rat2d, rat2di, tauold
write (4, * ) aveci, avecd, avec2d, avec2di
write (4, * ) avec0, sig0, v0, e_start, etot
write (4, * ) acu, ack, acp, acpv, avu, avk, avp, avpv, sigmamet
write (4, * ) istep, nzero, ntimes
close (unit = 4, status = 'keep')
WRITE( stdout,'(A3,3X,3F14.9)') ( atm(ityp(na)), tau(:,na), na = 1, nat )
WRITE( stdout, '(/5X,"Ekin = ",F14.8," Ryd T = ",F6.1," K ", &
& " Etot = ",F14.8)') ekint, tnew, edyn + e_start
!
! find alpha & beta for charge-density and wfc extrapolation (uptate_pot
! ... save MD history on file
!
if (calc(1:1).eq.'m') call find_alpha_and_beta(nat,tau,tauold,alpha0,beta0)
CALL seqopn( 4, TRIM( prefix )//'.md', 'FORMATTED', exst )
!
! Deallocate
WRITE(4,*) rati, ratd, rat2d, rat2di, tauold
WRITE(4,*) aveci, avecd, avec2d, avec2di
WRITE(4,*) avec0, sig0, v0, e_start, etot
WRITE(4,*) acu, ack, acp, acpv, avu, avk, avp, avpv, sigmamet
WRITE(4,*) istep, nzero, ntimes
!
deallocate (rat, rati, ratd, rat2d, rat2di, tauold)
CLOSE( UNIT = 4, STATUS = 'KEEP' )
!
return
101 format(1x,4d12.5,i6)
103 format(1x,6d12.5,i6)
104 format(1x,2d12.5,i6)
105 format(1x,3d12.5,i6)
end subroutine vcsmd
subroutine delete_if_present(filename)
USE io_global, ONLY : stdout
character (len=*) :: filename
logical :: exst, opnd
integer :: iunit
inquire (file = filename, exist = exst)
if (exst) then
do iunit = 99, 1, - 1
inquire (unit = iunit, opened = opnd)
if (.not.opnd) goto 10
enddo
call errore ('delete_if_present', 'free unit not found?!?', 1)
10 continue
open (unit=iunit, file= filename , status='old')
close (unit=iunit, status = 'delete')
WRITE( stdout,*) 'WARNING: ',filename,' file was present; old file deleted '
end if
end subroutine
! ... Deallocate
!
DEALLOCATE( rat, rati, ratd, rat2d, rat2di, tauold )
!
RETURN
!
101 FORMAT(1X,4D12.5,I6)
103 FORMAT(1X,6D12.5,I6)
104 FORMAT(1X,2D12.5,I6)
105 FORMAT(1X,3D12.5,I6)
!
END SUBROUTINE vcsmd