From eac4123a5b19158733c5343201d535d856d0418b Mon Sep 17 00:00:00 2001 From: sbraccia Date: Sun, 25 Jun 2006 23:15:20 +0000 Subject: [PATCH] Added to CP the possibility to remove the total torque acting on the system (useful to simulate an isolated system). Documentation updated. C.S. git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3205 c92efa57-630b-4861-b058-cf58834340f0 --- CPV/cpr.f90 | 31 ++++++++++++++++++------------- CPV/input.f90 | 37 ++++++++++++++++++------------------- Doc/INPUT_CP | 15 +++++++++++++++ 3 files changed, 51 insertions(+), 32 deletions(-) diff --git a/CPV/cpr.f90 b/CPV/cpr.f90 index 0fc9b5067..a43c702f8 100644 --- a/CPV/cpr.f90 +++ b/CPV/cpr.f90 @@ -13,14 +13,13 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out ) ! USE kinds, ONLY : DP USE constants, ONLY : bohr_radius_angs, amu_au - USE control_flags, ONLY : iprint, isave, thdyn, tpre, & - iprsta, tfor, tvlocw, & - taurdr, tprnfor, tsdc, lconstrain, lwf, & - lneb, lcoarsegrained, ndr, ndw, nomore, & - tsde, tortho, tnosee, tnosep, trane, & - tranp, tsdp, tcp, tcap, ampre, amprp, & - tnoseh, tolp, ortho_eps, ortho_max, & - printwfc + USE control_flags, ONLY : iprint, isave, thdyn, tpre, iprsta, & + tfor, tvlocw,remove_rigid_rot, taurdr, & + tprnfor, tsdc, lconstrain, lwf, lneb, & + lcoarsegrained, ndr, ndw, nomore, tsde, & + tortho, tnosee, tnosep, trane, tranp, & + tsdp, tcp, tcap, ampre, amprp, tnoseh, & + tolp, ortho_eps, ortho_max, printwfc USE core, ONLY : nlcc_any, rhoc USE uspp_param, ONLY : nhm, nh USE cvan, ONLY : nvb, ish @@ -156,7 +155,7 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out ) ! ! ... forces on ions ! - REAL(DP) :: maxfion + REAL(DP) :: maxfion, fion_tot(3) ! ! ... work variables ! @@ -168,11 +167,9 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out ) REAL(DP) :: delta_etot REAL(DP) :: ftmp, enb, enbi INTEGER :: is, nacc, ia, j, iter, i, isa, ipos - INTEGER :: k, ii, l, m, iss - ! + INTEGER :: k, ii, l, m, iss REAL(DP) :: hgamma(3,3), temphh(3,3) REAL(DP) :: fcell(3,3) - ! REAL(DP) :: stress_gpa(3,3), thstress(3,3) ! REAL(DP), ALLOCATABLE :: usrt_tau0(:,:), usrt_taup(:,:), usrt_fion(:,:) @@ -181,9 +178,10 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out ) REAL(DP), ALLOCATABLE :: tauw(:,:) ! temporary array used to printout positions CHARACTER(LEN=3) :: labelw( nat ) - ! for force_pairing + ! for force_pairing INTEGER :: nspin_sub ! + ! dt2bye = dt2 / emass etot_out = 0.D0 enow = 1.D9 @@ -335,6 +333,13 @@ SUBROUTINE cprmain( tau_out, fion_out, etot_out ) ! IF ( lwf ) CALL ef_force( fion, na, nsp, zv ) ! + fion_tot(:) = SUM( fion(:,:), DIM = 2 ) / DBLE( nat ) + ! + FORALL( ia = 1:nat ) fion(:,ia) = fion(:,ia) - fion_tot(:) + ! + IF ( remove_rigid_rot ) & + CALL remove_tot_torque( nat, tau0, pmass(ityp(ind_srt(:))), fion ) + ! IF ( lconstrain ) THEN ! IF ( ionode ) THEN diff --git a/CPV/input.f90 b/CPV/input.f90 index fedb30ec5..27009d157 100644 --- a/CPV/input.f90 +++ b/CPV/input.f90 @@ -261,6 +261,7 @@ MODULE input etot_maxiter_ => etot_maxiter, & forc_maxiter_ => forc_maxiter USE control_flags, ONLY : force_pairing_ => force_pairing + USE control_flags, ONLY : remove_rigid_rot_ => remove_rigid_rot ! ! ... Other modules ! @@ -278,8 +279,9 @@ MODULE input tefield2_ => tefield2, & epol2_ => epol2, & efield2_ => efield2 - - + ! + USE cvan, ONLY: nvb + ! USE input_parameters, ONLY: & electron_dynamics, electron_damping, diis_rot, electron_temperature, & ion_dynamics, ekin_conv_thr, etot_conv_thr, forc_conv_thr, ion_maxstep,& @@ -292,9 +294,7 @@ MODULE input tdipole_card, toptical_card, tnewnfi_card, newnfi_card, & ampre, nstep, restart_mode, ion_positions, startingwfc, printwfc, & orthogonalization, electron_velocities, nat, if_pos, phase_space, & - tefield, epol, efield, tefield2, epol2, efield2 - - USE cvan, ONLY: nvb + tefield, epol, efield, tefield2, epol2, efield2, remove_rigid_rot ! IMPLICIT NONE ! @@ -316,16 +316,15 @@ MODULE input etot_conv_thr_ = etot_conv_thr forc_conv_thr_ = forc_conv_thr ekin_maxiter_ = electron_maxstep - - - tefield_ = tefield - epol_ = epol - efield_ = efield - - tefield2_ = tefield2 - epol2_ = epol2 - efield2_ = efield2 - + ! + remove_rigid_rot_ = remove_rigid_rot + ! + tefield_ = tefield + epol_ = epol + efield_ = efield + tefield2_ = tefield2 + epol2_ = epol2 + efield2_ = efield2 ! ! ... Set internal time step variables ( delt, twodelt, dt2 ... ) ! @@ -354,8 +353,8 @@ MODULE input trhor_ = ( TRIM( calculation ) == 'nscf' ) tvlocw_ = ( TRIM( disk_io ) == 'high' ) ! warning this is not working ! for now use reduce_io in CP to specify if the charge density is saved -! reduce_io_ = ( TRIM( disk_io ) == 'low' ) - reduce_io_ = .not.( TRIM( disk_io ) == 'high' ) + ! + reduce_io_ = .NOT.( TRIM( disk_io ) == 'high' ) ! SELECT CASE( TRIM( verbosity ) ) CASE( 'minimal' ) @@ -661,6 +660,8 @@ MODULE input CALL errore(' control_flags ',' unknown ion_dynamics '//TRIM(ion_dynamics), 1 ) END SELECT + + IF ( ANY( if_pos(:,1:nat) == 0 ) ) lfixatom = .TRUE. ! ... Ionic Temperature @@ -806,8 +807,6 @@ MODULE input force_pairing_ = force_pairing ! - ! - ! IF( ( nvb > 0 ) .and. ( program_name == 'FPMD' ) ) & CALL errore(' iosys ',' USPP not yet implemented in FPMD ',1) diff --git a/Doc/INPUT_CP b/Doc/INPUT_CP index ae901f9fa..94bc588cf 100644 --- a/Doc/INPUT_CP +++ b/Doc/INPUT_CP @@ -564,6 +564,21 @@ amprp(i) REAL ( default = 0.D0 ) greasp REAL ( default = 1.D0 ) same as "grease", for ionic damped dynamics. +remove_rigid_rot + LOGICAL ( default = .FALSE. ) + this keyword is useful when simulating the dynamics and/or the + thermodynamics of an isolated system. If set to true the total + torque of the internal forces is set to zero by adding new forces + that compensate the spurious interaction with the periodic + images. This allowes for the use of smaller supercells. + BEWARE: since the potential energy is no longer consistent with + the forces (it still contains the spurious interaction with the + repeated images), the total energy is not conserved anymore. + However the dynamical and thermodynamical properties should be + in closer agreement with those of an isolated system. + Also the final energy of a structural relaxation will be higher, + but the relaxation itself should be faster. + ! ! ... keywords used only in NEB calculations !