From 0224c2a9ad059196b64bb4881e8267cc350566a5 Mon Sep 17 00:00:00 2001 From: sbraccia Date: Thu, 26 Feb 2004 11:50:36 +0000 Subject: [PATCH] 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 --- Modules/control_flags.f90 | 4 +- Modules/input_parameters.f90 | 8 +- Modules/parser.f90 | 118 +++-- Modules/read_cards.f90 | 2 +- Modules/read_namelists.f90 | 2 +- PW/Makefile | 5 +- PW/add_vuspsi.f90 | 75 ++- PW/clean_pw.f90 | 6 + PW/dynamics.f90 | 877 +++++++++++++++++------------------ PW/input.f90 | 157 ++++--- PW/move_ions.f90 | 494 ++++++++++++-------- PW/update_pot.f90 | 19 +- PW/vcsmd.f90 | 871 ++++++++++++++-------------------- 13 files changed, 1346 insertions(+), 1292 deletions(-) diff --git a/Modules/control_flags.f90 b/Modules/control_flags.f90 index ca719509e..660b8367a 100644 --- a/Modules/control_flags.f90 +++ b/Modules/control_flags.f90 @@ -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 diff --git a/Modules/input_parameters.f90 b/Modules/input_parameters.f90 index b49532095..a43627778 100644 --- a/Modules/input_parameters.f90 +++ b/Modules/input_parameters.f90 @@ -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 diff --git a/Modules/parser.f90 b/Modules/parser.f90 index 5ccc713d0..f4e9af027 100644 --- a/Modules/parser.f90 +++ b/Modules/parser.f90 @@ -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) !----------------------------------------------------------------------- diff --git a/Modules/read_cards.f90 b/Modules/read_cards.f90 index 46b728375..137dfa111 100644 --- a/Modules/read_cards.f90 +++ b/Modules/read_cards.f90 @@ -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 diff --git a/Modules/read_namelists.f90 b/Modules/read_namelists.f90 index 3934cce77..8ba17e871 100644 --- a/Modules/read_namelists.f90 +++ b/Modules/read_namelists.f90 @@ -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 diff --git a/PW/Makefile b/PW/Makefile index 50e2db624..21c1e79c9 100644 --- a/PW/Makefile +++ b/PW/Makefile @@ -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 \ diff --git a/PW/add_vuspsi.f90 b/PW/add_vuspsi.f90 index 5c4718b0a..b7b85aa84 100644 --- a/PW/add_vuspsi.f90 +++ b/PW/add_vuspsi.f90 @@ -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 diff --git a/PW/clean_pw.f90 b/PW/clean_pw.f90 index ee33e2e11..198e0aa5c 100644 --- a/PW/clean_pw.f90 +++ b/PW/clean_pw.f90 @@ -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 diff --git a/PW/dynamics.f90 b/PW/dynamics.f90 index 6d9dc3130..c2fda368b 100644 --- a/PW/dynamics.f90 +++ b/PW/dynamics.f90 @@ -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 diff --git a/PW/input.f90 b/PW/input.f90 index 12b8ac06d..ad10a65c5 100644 --- a/PW/input.f90 +++ b/PW/input.f90 @@ -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 diff --git a/PW/move_ions.f90 b/PW/move_ions.f90 index 6523f7f06..826c512e0 100644 --- a/PW/move_ions.f90 +++ b/PW/move_ions.f90 @@ -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 diff --git a/PW/update_pot.f90 b/PW/update_pot.f90 index 86410d655..f4c3d4c9f 100644 --- a/PW/update_pot.f90 +++ b/PW/update_pot.f90 @@ -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 ) diff --git a/PW/vcsmd.f90 b/PW/vcsmd.f90 index 5f9ec4c69..164ede69e 100644 --- a/PW/vcsmd.f90 +++ b/PW/vcsmd.f90 @@ -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