1! 2! Copyright (C) 2001-2011 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9!---------------------------------------------------------------------------- 10SUBROUTINE forces() 11 !---------------------------------------------------------------------------- 12 ! 13 ! ... This routine is a driver routine which computes the forces 14 ! ... acting on the atoms. The complete expression of the forces 15 ! ... contains four parts which are computed by different routines: 16 ! 17 ! ... a) force_lc, local contribution to the forces 18 ! ... b) force_cc, contribution due to NLCC 19 ! ... c) force_ew, contribution due to the electrostatic ewald term 20 ! ... d) force_us, contribution due to the non-local potential 21 ! ... e) force_corr, correction term for incomplete self-consistency 22 ! ... f) force_hub, contribution due to the Hubbard term 23 ! ... g) force_london, semi-empirical correction for dispersion forces 24 ! 25 ! 26 USE kinds, ONLY : DP 27 USE io_global, ONLY : stdout 28 USE cell_base, ONLY : at, bg, alat, omega 29 USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau, zv, amass, extfor 30 USE fft_base, ONLY : dfftp 31 USE gvect, ONLY : ngm, gstart, ngl, nl, igtongl, g, gg, gcutm 32 USE lsda_mod, ONLY : nspin 33 USE symme, ONLY : symvector 34 USE vlocal, ONLY : strf, vloc 35 USE force_mod, ONLY : force, lforce 36 USE scf, ONLY : rho 37 USE ions_base, ONLY : if_pos 38 USE ldaU, ONLY : lda_plus_u, U_projection 39 USE extfield, ONLY : tefield, forcefield 40 USE control_flags, ONLY : gamma_only, remove_rigid_rot, textfor, & 41 iverbosity, llondon 42#ifdef __ENVIRON 43 USE environ_base, ONLY : do_environ, env_static_permittivity, rhopol 44 USE fft_interfaces, ONLY : fwfft 45#endif 46 USE bp, ONLY : lelfield, gdir, l3dstring, efield_cart, & 47 efield_cry,efield 48 USE uspp, ONLY : okvan 49 USE martyna_tuckerman, ONLY: do_comp_mt, wg_corr_force 50 USE london_module, ONLY : force_london 51 ! 52 IMPLICIT NONE 53 ! 54 REAL(DP), ALLOCATABLE :: forcenl(:,:), & 55 forcelc(:,:), & 56 forcecc(:,:), & 57 forceion(:,:), & 58 force_disp(:,:),& 59 force_mt(:,:), & 60 forcescc(:,:), & 61 forces_bp_efield(:,:), & 62 forceh(:,:) 63 ! nonlocal, local, core-correction, ewald, scf correction terms, and hubbard 64#ifdef __ENVIRON 65 REAL(DP), ALLOCATABLE :: force_environ(:,:) 66 COMPLEX(DP), ALLOCATABLE :: aux(:) 67#endif 68 69 REAL(DP) :: sumfor, sumscf, sum_mm 70 REAL(DP),PARAMETER :: eps = 1.e-12_dp 71 INTEGER :: ipol, na 72 ! counter on polarization 73 ! counter on atoms 74 ! 75 ! 76 CALL start_clock( 'forces' ) 77 ! 78 ALLOCATE( forcenl( 3, nat ), forcelc( 3, nat ), forcecc( 3, nat ), & 79 forceh( 3, nat ), forceion( 3, nat ), forcescc( 3, nat ) ) 80 ! 81 forcescc(:,:) = 0.D0 82 forceh(:,:) = 0.D0 83 ! 84 WRITE( stdout, '(/,5x,"Forces acting on atoms (Ry/au):", / )') 85 ! 86 ! ... The nonlocal contribution is computed here 87 ! 88 CALL force_us( forcenl ) 89 ! 90 ! ... The local contribution 91 ! 92 CALL force_lc( nat, tau, ityp, alat, omega, ngm, ngl, igtongl, & 93 g, rho%of_r, nl, nspin, gstart, gamma_only, vloc, & 94 forcelc ) 95 ! 96 ! ... The NLCC contribution 97 ! 98 CALL force_cc( forcecc ) 99 ! 100 ! ... The Hubbard contribution 101 ! (included by force_us if using beta as local projectors) 102 ! 103 IF ( lda_plus_u .AND. U_projection.NE.'pseudo' ) CALL force_hub( forceh ) 104 ! 105 ! ... The ionic contribution is computed here 106 ! 107 CALL force_ew( alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, & 108 gg, ngm, gstart, gamma_only, gcutm, strf, forceion ) 109 ! 110 ! ... the semi-empirical dispersion correction 111 ! 112 IF ( llondon ) THEN 113 ! 114 ALLOCATE ( force_disp ( 3 , nat ) ) 115 force_disp ( : , : ) = 0.0_DP 116 force_disp = force_london( alat , nat , ityp , at , bg , tau ) 117 ! 118 END IF 119 120 ! 121 ! ... The SCF contribution 122 ! 123 CALL force_corr( forcescc ) 124 ! 125 IF (do_comp_mt) THEN 126 ! 127 ALLOCATE ( force_mt ( 3 , nat ) ) 128#ifdef __ENVIRON 129 IF ( do_environ .AND. env_static_permittivity .GT. 1.D0 ) THEN 130 ALLOCATE( aux( dfftp%nnr ) ) 131 aux(:) = CMPLX(rhopol( : ),0.D0,kind=dp) 132 CALL fwfft ('Dense', aux, dfftp) 133 rho%of_g(:,1) = rho%of_g(:,1) + aux(nl(:)) 134 ENDIF 135#endif 136 CALL wg_corr_force( omega, nat, ntyp, ityp, ngm, g, tau, zv, strf, & 137 nspin, rho%of_g, force_mt ) 138#ifdef __ENVIRON 139 IF ( do_environ .AND. env_static_permittivity .GT. 1.D0 ) THEN 140 rho%of_g(:,1) = rho%of_g(:,1) - aux(nl(:)) 141 DEALLOCATE(aux) 142 ENDIF 143#endif 144 145 END IF 146#ifdef __ENVIRON 147 IF (do_environ) THEN 148 ! 149 ! ... The external environment contribution 150 ! 151 ALLOCATE ( force_environ ( 3 , nat ) ) 152 force_environ = 0.0_DP 153 ! 154 ! ... Computes here the solvent contribution 155 ! 156 IF ( env_static_permittivity .GT. 1.D0 ) & 157 CALL force_lc( nat, tau, ityp, alat, omega, ngm, ngl, igtongl, & 158 g, rhopol, nl, 1, gstart, gamma_only, vloc, & 159 force_environ ) 160 ! 161 ! ... Add the other environment contributions 162 ! 163 CALL calc_fenviron( dfftp%nnr, nat, force_environ ) 164 ! 165 END IF 166 ! 167#endif 168 ! 169 ! Berry's phase electric field terms 170 ! 171 if(lelfield) then 172 ALLOCATE ( forces_bp_efield (3,nat) ) 173 forces_bp_efield(:,:)=0.d0 174 if(.not.l3dstring) then 175 if(okvan) call forces_us_efield(forces_bp_efield,gdir,efield) 176 call forces_ion_efield(forces_bp_efield,gdir,efield) 177 else 178 if(okvan)then 179 do ipol=1,3 180 call forces_us_efield(forces_bp_efield,ipol,efield_cry(ipol)) 181 enddo 182 endif 183 do ipol=1,3 184 call forces_ion_efield(forces_bp_efield,ipol,efield_cart(ipol)) 185 enddo 186 endif 187 endif 188 ! 189 ! ... here we sum all the contributions and compute the total force acting 190 ! ... on the crstal 191 ! 192 DO ipol = 1, 3 193 ! 194 sumfor = 0.D0 195 ! 196 DO na = 1, nat 197 ! 198 force(ipol,na) = forcenl(ipol,na) + & 199 forceion(ipol,na) + & 200 forcelc(ipol,na) + & 201 forcecc(ipol,na) + & 202 forceh(ipol,na) + & 203 forcescc(ipol,na) 204 ! 205 IF ( llondon ) force(ipol,na) = force(ipol,na) + force_disp(ipol,na) 206 IF ( tefield ) force(ipol,na) = force(ipol,na) + forcefield(ipol,na) 207 IF (lelfield) force(ipol,na) = force(ipol,na) + forces_bp_efield(ipol,na) 208 IF (do_comp_mt)force(ipol,na) = force(ipol,na) + force_mt(ipol,na) 209! DCC 210! IF (do_comp) force(ipol,na) = force(ipol,na) + force_vcorr(ipol,na) 211#ifdef __ENVIRON 212 IF (do_environ) force(ipol,na) = force(ipol,na) + force_environ(ipol,na) 213#endif 214 215 sumfor = sumfor + force(ipol,na) 216 ! 217 END DO 218 ! 219 ! ... impose total force = 0 220 ! 221 DO na = 1, nat 222 ! 223 force(ipol,na) = force(ipol,na) - sumfor / DBLE( nat ) 224 ! 225 END DO 226 ! 227#ifdef __MS2 228 ! 229 ! ... impose total force of the quantum subsystem /= 0 230 ! 231 DO na = 1, nat 232 ! 233 force(ipol,na) = force(ipol,na) + sumfor / DBLE( nat ) 234 ! 235 END DO 236 ! 237#endif 238 ! 239 END DO 240 ! 241 ! ... resymmetrize (should not be needed, but ...) 242 ! 243 CALL symvector ( nat, force ) 244 ! 245 IF ( remove_rigid_rot ) & 246 CALL remove_tot_torque( nat, tau, amass(ityp(:)), force ) 247 ! 248 IF( textfor ) force(:,:) = force(:,:) + extfor(:,:) 249 ! 250 ! ... call void routine for user define/ plugin patches on forces 251 ! 252 ! plugin should be called after stress has been computed 253 ! Look in pwscf.f90 254 ! CALL plugin_forces() 255 ! 256 ! ... write on output the forces 257 ! 258 DO na = 1, nat 259 ! 260 WRITE( stdout, 9035) na, ityp(na), force(:,na) 261 ! 262 END DO 263 ! 264 ! ... forces on fixed coordinates are set to zero ( C.S. 15/10/2003 ) 265 ! 266 force(:,:) = force(:,:) * DBLE( if_pos ) 267 forcescc(:,:) = forcescc(:,:) * DBLE( if_pos ) 268 ! 269 IF ( iverbosity > 0 ) THEN 270 IF ( do_comp_mt ) THEN 271 WRITE( stdout, '(5x,"The Martyna-Tuckerman correction term to forces")') 272 DO na = 1, nat 273 WRITE( stdout, 9035) na, ityp(na), ( force_mt(ipol,na), ipol = 1, 3 ) 274 END DO 275 END IF 276 ! 277 WRITE( stdout, '(5x,"The non-local contrib. to forces")') 278 DO na = 1, nat 279 WRITE( stdout, 9035) na, ityp(na), ( forcenl(ipol,na), ipol = 1, 3 ) 280 END DO 281 WRITE( stdout, '(5x,"The ionic contribution to forces")') 282 DO na = 1, nat 283 WRITE( stdout, 9035) na, ityp(na), ( forceion(ipol,na), ipol = 1, 3 ) 284 END DO 285 WRITE( stdout, '(5x,"The local contribution to forces")') 286 DO na = 1, nat 287 WRITE( stdout, 9035) na, ityp(na), ( forcelc(ipol,na), ipol = 1, 3 ) 288 END DO 289 WRITE( stdout, '(5x,"The core correction contribution to forces")') 290 DO na = 1, nat 291 WRITE( stdout, 9035) na, ityp(na), ( forcecc(ipol,na), ipol = 1, 3 ) 292 END DO 293 WRITE( stdout, '(5x,"The Hubbard contrib. to forces")') 294 DO na = 1, nat 295 WRITE( stdout, 9035) na, ityp(na), ( forceh(ipol,na), ipol = 1, 3 ) 296 END DO 297 WRITE( stdout, '(5x,"The SCF correction term to forces")') 298 DO na = 1, nat 299 WRITE( stdout, 9035) na, ityp(na), ( forcescc(ipol,na), ipol = 1, 3 ) 300 END DO 301#ifdef __ENVIRON 302 IF ( do_environ ) THEN 303 WRITE( stdout, '(5x,"The external environment correction to forces")') 304 DO na = 1, nat 305 WRITE( stdout, 9035) na, ityp(na), ( force_environ(ipol,na), ipol = 1, 3 ) 306 END DO 307 END IF 308#endif 309 ! 310 IF ( llondon) THEN 311 WRITE( stdout, '(/,5x,"Dispersion contribution to forces:")') 312 DO na = 1, nat 313 WRITE( stdout, 9035) na, ityp(na), (force_disp(ipol,na), ipol = 1, 3) 314 END DO 315 END IF 316 ! 317 END IF 318 ! 319 sumfor = 0.D0 320 sumscf = 0.D0 321 ! 322 DO na = 1, nat 323 ! 324 sumfor = sumfor + force(1,na)**2 + force(2,na)**2 + force(3,na)**2 325 sumscf = sumscf + forcescc(1,na)**2 + forcescc(2,na)**2+ forcescc(3,na)**2 326 ! 327 END DO 328 ! 329 sumfor = SQRT( sumfor ) 330 sumscf = SQRT( sumscf ) 331 ! 332 WRITE( stdout, '(/5x,"Total force = ",F12.6,5X, & 333 & "Total SCF correction = ",F12.6)') sumfor, sumscf 334 ! 335 IF ( llondon .AND. iverbosity > 0 ) THEN 336 ! 337 sum_mm = 0.D0 338 DO na = 1, nat 339 sum_mm = sum_mm + & 340 force_disp(1,na)**2 + force_disp(2,na)**2 + force_disp(3,na)**2 341 END DO 342 sum_mm = SQRT( sum_mm ) 343 WRITE ( stdout, '(/,5x, "Total Dispersion Force = ",F12.6)') sum_mm 344 ! 345 END IF 346 ! 347 DEALLOCATE( forcenl, forcelc, forcecc, forceh, forceion, forcescc ) 348 IF ( llondon ) DEALLOCATE ( force_disp ) 349 IF ( lelfield ) DEALLOCATE ( forces_bp_efield ) 350 ! 351 lforce = .TRUE. 352 ! 353 CALL stop_clock( 'forces' ) 354 ! 355 IF ( ( sumfor < 10.D0*sumscf ) .AND. ( sumfor > eps ) ) & 356 WRITE( stdout,'(5x,"SCF correction compared to forces is large: ", & 357 & "reduce conv_thr to get better values")') 358 ! 359 IF(ALLOCATED(force_mt)) DEALLOCATE( force_mt ) 360#ifdef __ENVIRON 361 IF(ALLOCATED(force_environ)) DEALLOCATE( force_environ ) 362#endif 363 364 RETURN 365 ! 3669035 FORMAT(5X,'atom ',I4,' type ',I2,' force = ',3F14.8) 367 ! 368END SUBROUTINE forces 369