1! 2! Copyright (C) 2002-2010 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!=----------------------------------------------------------------------=! 10! 11! CP90 / FPMD common init subroutine 12! 13!=----------------------------------------------------------------------=! 14 15 subroutine init_dimensions( ) 16 17 ! 18 ! initialize G-vectors and related quantities 19 ! 20 21 USE kinds, ONLY : dp 22 USE constants, ONLY : tpi 23 use io_global, only : stdout, ionode 24 use control_flags, only : gamma_only, iverbosity 25 use cell_base, only : ainv, at, omega, alat 26 use small_box, only : small_box_set 27 use smallbox_grid_dim, only : smallbox_grid_init,smallbox_grid_info 28 USE fft_types, ONLY : fft_type_init 29 use ions_base, only : nat 30 USE recvec_subs, ONLY : ggen, ggens 31 USE gvect, ONLY : mill_g, eigts1,eigts2,eigts3, g, gg, & 32 ecutrho, gcutm, gvect_init, mill, & 33 ig_l2g, gstart, ngm, ngm_g, gshells 34 use gvecs, only : gcutms, gvecs_init, ngms 35 use gvecw, only : gkcut, gvecw_init, g2kin_init 36 USE smallbox_subs, ONLY : ggenb 37 USE fft_base, ONLY : dfftp, dffts, dfftb, fft_base_info 38 USE fft_smallbox, ONLY : cft_b_omp_init 39 USE fft_base, ONLY : smap 40 USE control_flags, ONLY : gamma_only, smallmem 41 USE electrons_module, ONLY : bmeshset 42 USE electrons_base, ONLY : distribute_bands 43 USE problem_size, ONLY : cpsizes 44 USE mp_bands, ONLY : me_bgrp, root_bgrp, nproc_bgrp, nbgrp, & 45 my_bgrp_id, intra_bgrp_comm, ntask_groups 46 USE uspp, ONLY : okvan, nlcc_any 47 USE input_parameters, ONLY : ref_cell, ref_alat 48 use cell_base, ONLY : ref_at, ref_bg 49 USE exx_module, ONLY : h_init 50 USE command_line_options, ONLY : nmany_ 51 52 implicit none 53! 54 integer :: i 55 real(dp) :: rat1, rat2, rat3 56 real(dp) :: bg(3,3), tpiba2 57 integer :: ng_, ngs_, ngm_ , ngw_, nyfft_ 58#if defined(__MPI) 59 LOGICAL :: lpara = .true. 60#else 61 LOGICAL :: lpara = .false. 62#endif 63 64 CALL start_clock( 'init_dim' ) 65 66 tpiba2 = ( tpi / alat ) ** 2 67 IF( ionode ) THEN 68 WRITE( stdout, 100 ) 69 100 FORMAT( //, & 70 3X,'Simulation dimensions initialization',/, & 71 3X,'------------------------------------' ) 72 END IF 73 ! 74 ! ... Initialize bands indexes for parallel linear algebra 75 ! ... (distribute bands to processors) 76 ! 77 CALL bmeshset( ) 78 ! 79 ! ... cell dimensions and lattice vectors 80 ! ... note that at are in alat units 81 82 call recips( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) ) 83 84 ! bg(:,1), bg(:,2), bg(:,3) are the basis vectors, in 85 ! 2pi/alat units, generating the reciprocal lattice 86 87 ! Store the cell parameter from the input file. Used in exx_module ... 88 h_init=at*alat 89 90 ! ... Initialize FFT real-space grids and small box grid 91 nyfft_ = ntask_groups 92 dffts%has_task_groups = (ntask_groups >1) 93 dfftp%has_task_groups = .FALSE. 94 lpara = ( nproc_bgrp > 1 ) 95 ! 96 IF ( ref_cell ) THEN 97 ! 98 CALL recips( ref_at(1,1), ref_at(1,2), ref_at(1,3), ref_bg(1,1), ref_bg(1,2), ref_bg(1,3) ) 99 ! 100 WRITE( stdout,'(3X,"Reference Cell is Used to Initialize FFT Real-space Grids")' ) 101 WRITE( stdout,'(3X,"Reference Cell alat =",F14.8,1X,"A.U.")' ) ref_alat 102 WRITE( stdout,'(3X,"ref_cell_a1 =",1X,3f14.8,3x,"ref_cell_b1 =",3f14.8)') ref_at(:,1)*ref_alat,ref_bg(:,1)/ref_alat 103 WRITE( stdout,'(3X,"ref_cell_a2 =",1X,3f14.8,3x,"ref_cell_b2 =",3f14.8)') ref_at(:,2)*ref_alat,ref_bg(:,2)/ref_alat 104 WRITE( stdout,'(3X,"ref_cell_a3 =",1X,3f14.8,3x,"ref_cell_b3 =",3f14.8)') ref_at(:,3)*ref_alat,ref_bg(:,3)/ref_alat 105 ! 106 CALL fft_type_init( dffts, smap, "wave", gamma_only, lpara, intra_bgrp_comm, ref_at, ref_bg, & 107 gkcut, nyfft=nyfft_, nmany=nmany_ ) 108 CALL fft_type_init( dfftp, smap, "rho", gamma_only, lpara, intra_bgrp_comm, ref_at, ref_bg, & 109 gcutm, nyfft=nyfft_, nmany=nmany_ ) 110 ! 111 ELSE 112 ! 113 CALL fft_type_init( dffts, smap, "wave", gamma_only, lpara, intra_bgrp_comm, at, bg, & 114 gkcut, nyfft=nyfft_, nmany=nmany_ ) 115 CALL fft_type_init( dfftp, smap, "rho", gamma_only, lpara, intra_bgrp_comm, at, bg, & 116 gcutm, nyfft=nyfft_, nmany=nmany_ ) 117 ! 118 END IF 119 ! define the clock labels ( this enables the corresponding fft too ! ) 120 dffts%rho_clock_label = 'ffts' ; dffts%wave_clock_label = 'fftw' 121 dfftp%rho_clock_label = 'fft' 122 ! 123 ! 124 CALL smallbox_grid_init( dfftp, dfftb ) 125 126 IF( ionode ) THEN 127 128 WRITE( stdout,210) 129210 format(/,3X,'unit vectors of full simulation cell',& 130 &/,3X,'in real space:',25x,'in reciprocal space (units 2pi/alat):') 131 WRITE( stdout,'(3X,I1,1X,3f10.4,10x,3f10.4)') 1,at(:,1)*alat,bg(:,1) 132 WRITE( stdout,'(3X,I1,1X,3f10.4,10x,3f10.4)') 2,at(:,2)*alat,bg(:,2) 133 WRITE( stdout,'(3X,I1,1X,3f10.4,10x,3f10.4)') 3,at(:,3)*alat,bg(:,3) 134 135 END IF 136 ! 137 do i=1,3 138 ainv(1,i)=bg(i,1)/alat 139 ainv(2,i)=bg(i,2)/alat 140 ainv(3,i)=bg(i,3)/alat 141 end do 142 143 ! 144 ! ainv is transformation matrix from cartesian to crystal coordinates 145 ! if r=x1*a1+x2*a2+x3*a3 => x(i)=sum_j ainv(i,j)r(j) 146 ! Note that ainv is really the inverse of a=(a1,a2,a3) 147 ! (but only if the axis triplet is right-handed, otherwise 148 ! for a left-handed triplet, ainv is minus the inverse of a) 149 ! 150 CALL fft_base_info( ionode, stdout ) 151 ngw_ = dffts%ngw 152 ngs_ = dffts%ngm 153 ngm_ = dfftp%ngm 154 155 ! 156 ! ... Initialize reciprocal space local and global dimensions 157 ! NOTE in a parallel run ngm_ , ngw_ , ngs_ here are the 158 ! local number of reciprocal vectors 159 ! 160 CALL gvect_init ( ngm_ , intra_bgrp_comm ) 161 CALL gvecs_init ( ngs_ , intra_bgrp_comm ) 162 ! 163 ! ... Print real-space grid dimensions 164 ! 165 CALL realspace_grids_info ( dfftp, dffts ) 166 CALL smallbox_grid_info ( dfftb ) 167 ! 168 ! ... generate g-space vectors (dense and smooth grid) 169 ! ... call to gshells generates gl, igtongl used in vdW-DF functional 170 ! 171 IF ( ref_cell ) THEN 172 ! 173 WRITE( stdout,'(/,3X,"Reference Cell is Used to Initialize Reciprocal Space Mesh")' ) 174 WRITE( stdout,'(3X,"Reference Cell alat =",F14.8,1X,"A.U.")' ) ref_alat 175 ! 176 IF( smallmem ) THEN 177 CALL ggen( dfftp, gamma_only, ref_at, ref_bg, gcutm, ngm_g, ngm, & 178 g, gg, mill, ig_l2g, gstart, no_global_sort = .TRUE. ) 179 ELSE 180 CALL ggen( dfftp, gamma_only, ref_at, ref_bg, gcutm, ngm_g, ngm, & 181 g, gg, mill, ig_l2g, gstart ) 182 END IF 183 CALL ggens( dffts, gamma_only, ref_at, g, gg, mill, gcutms, ngms ) 184 ! 185 ELSE 186 ! 187 IF( smallmem ) THEN 188 CALL ggen( dfftp, gamma_only, at, bg, gcutm, ngm_g, ngm, & 189 g, gg, mill, ig_l2g, gstart, no_global_sort = .TRUE. ) 190 ELSE 191 CALL ggen( dfftp, gamma_only, at, bg, gcutm, ngm_g, ngm, & 192 g, gg, mill, ig_l2g, gstart ) 193 END IF 194 CALL ggens( dffts, gamma_only, at, g, gg, mill, gcutms, ngms ) 195 ! 196 END IF 197 198 CALL gshells (.TRUE.) 199 ! 200 ! ... allocate and generate (modified) kinetic energy 201 ! 202 CALL gvecw_init ( ngw_ , intra_bgrp_comm ) 203 CALL g2kin_init ( gg, tpiba2 ) 204 ! 205 ! global arrays are no more needed 206 ! 207 if( allocated( mill_g ) ) deallocate( mill_g ) 208 ! 209 ! allocate spaces for phases e^{-iG*tau_s} 210 ! 211 allocate( eigts1(-dfftp%nr1:dfftp%nr1,nat) ) 212 allocate( eigts2(-dfftp%nr2:dfftp%nr2,nat) ) 213 allocate( eigts3(-dfftp%nr3:dfftp%nr3,nat) ) 214 ! 215 ! small boxes 216 ! 217 IF ( dfftb%nr1 > 0 .AND. dfftb%nr2 > 0 .AND. dfftb%nr3 > 0 ) THEN 218 219 ! set the small box parameters 220 221 rat1 = DBLE( dfftb%nr1 ) / DBLE( dfftp%nr1 ) 222 rat2 = DBLE( dfftb%nr2 ) / DBLE( dfftp%nr2 ) 223 rat3 = DBLE( dfftb%nr3 ) / DBLE( dfftp%nr3 ) 224 ! 225 CALL small_box_set( alat, omega, at, rat1, rat2, rat3, tprint = .TRUE. ) 226 ! 227 ! generate small-box G-vectors, initialize FFT tables 228 ! 229 CALL ggenb ( ecutrho, iverbosity ) 230 ! 231#if defined _OPENMP 232 CALL cft_b_omp_init( dfftb%nr1, dfftb%nr2, dfftb%nr3 ) 233#endif 234 ELSE IF( okvan .OR. nlcc_any ) THEN 235 236 CALL errore( ' init_dimensions ', ' nr1b, nr2b, nr3b must be given for ultrasoft and core corrected pp ', 1 ) 237 238 END IF 239 240 ! ... distribute bands 241 242 CALL distribute_bands( nbgrp, my_bgrp_id ) 243 244 ! ... printout g vector distribution summary 245 ! 246 CALL gmeshinfo() 247 ! 248 ! CALL cpsizes( ) Maybe useful 249 ! 250 ! Flush stdout 251 ! 252 FLUSH( stdout ) 253 ! 254 CALL stop_clock( 'init_dim' ) 255 ! 256 return 257 258 CONTAINS 259 260 SUBROUTINE fft_extra_info() 261 INTEGER :: ir 262 IF(ionode) THEN 263 WRITE(stdout,*) 'FFT parallelization for potentials' 264 WRITE(stdout,*) dfftp%nproc, dfftp%nproc2, dfftp%nproc3 265 WRITE(stdout,*) 'FFT parallelization for smooth grid' 266 WRITE(stdout,*) dffts%nproc, dffts%nproc2, dffts%nproc3 267 END IF 268 WRITE(1000+dfftp%mype,*) dfftp%nr1x, dfftp%nr2x 269 DO ir = 1, dfftp%nr1x * dfftp%nr2x 270 WRITE(1000+dfftp%mype,*) dfftp%isind( ir ) 271 END DO 272 END SUBROUTINE fft_extra_info 273 274 END SUBROUTINE init_dimensions 275 276!----------------------------------------------------------------------- 277 subroutine init_geometry ( ) 278!----------------------------------------------------------------------- 279! 280 USE kinds, ONLY: DP 281 use control_flags, only: iprint, thdyn, ndr, nbeg, tbeg 282 use io_global, only: stdout, ionode 283 use mp_global, only: nproc_bgrp, me_bgrp, intra_bgrp_comm, root_bgrp 284 use ions_base, only: na, nsp, nat, tau, if_pos 285 use cell_base, only: at, alat, r_to_s, cell_init, deth 286 use cell_base, only: ibrav, ainv, h, hold, tcell_base_init 287 USE ions_positions, ONLY: allocate_ions_positions, tau0, taus 288 use cp_restart_new, only: cp_read_cell 289 USE fft_base, ONLY: dfftb 290 USE fft_smallbox_type, ONLY: fft_box_allocate 291 USE cp_main_variables,ONLY: ht0, htm, taub 292 USE cp_interfaces, ONLY: newinit 293 USE constants, ONLY: amu_au 294 USE matrix_inversion 295 296 implicit none 297 ! 298 ! local 299 ! 300 integer :: i, j 301 real(DP) :: gvel(3,3), ht(3,3) 302 real(DP) :: xnhh0(3,3), xnhhm(3,3), vnhh(3,3), velh(3,3) 303 REAL(DP), ALLOCATABLE :: pmass(:), taus_srt( :, : ) 304 305 IF( .NOT. tcell_base_init ) & 306 CALL errore( ' init_geometry ', ' cell_base_init has not been call yet! ', 1 ) 307 308 IF( ionode ) THEN 309 WRITE( stdout, 100 ) 310 100 FORMAT( //, & 311 3X,'System geometry initialization',/, & 312 3X,'------------------------------' ) 313 END IF 314 315 ! Set ht0 and htm, cell at time t and t-dt 316 ! 317 CALL cell_init( alat, at, ht0 ) 318 CALL cell_init( alat, at, htm ) 319 320 CALL allocate_ions_positions( nsp, nat ) 321 ! 322 ! tau0 = initial positions, sorted wrt order read from input 323 ! taus = initial positions, scaled with the cell read from input 324 ! 325 tau0(:,:) = tau(:,:) 326 CALL r_to_s( tau, taus, nat, ainv ) 327 ! 328 ! Allocate box descriptor 329 ! 330 ALLOCATE( taub( 3, nat ) ) 331 ! 332 CALL fft_box_allocate( dfftb, me_bgrp, root_bgrp, nproc_bgrp, intra_bgrp_comm, nat ) 333 ! 334 ! if tbeg = .true. the geometry is given in the standard input even if 335 ! we are restarting a previous run 336 ! 337 if( ( nbeg > -1 ) .and. ( .not. tbeg ) ) then 338 ! 339 ! read only h and hold from restart file "ndr" 340 ! 341 CALL cp_read_cell( ndr, .TRUE., ht, hold, velh, gvel, xnhh0, xnhhm, vnhh ) 342 343 CALL cell_init( 't', ht0, ht ) 344 CALL cell_init( 't', htm, hold ) 345 ht0%hvel = velh ! set cell velocity 346 ht0%gvel = gvel 347 348 h = TRANSPOSE( ht ) 349 ht = TRANSPOSE( hold ) 350 hold = ht 351 ht = TRANSPOSE( velh ) 352 velh = ht 353 354 ! BS ... additional printing hold 355 WRITE(stdout, '(3X,"cell parameters read from restart file")') 356 WRITE( stdout,344) ibrav 357 WRITE(stdout, '(/,3X,"cell at current step : h(t)")') 358 do i=1,3 359 WRITE( stdout,345) (h(i,j),j=1,3) 360 enddo 361 WRITE(stdout, '(/,3X,"cell at previous step : h(t-dt)")') 362 do i=1,3 363 WRITE( stdout,345) (hold(i,j),j=1,3) 364 enddo 365 WRITE( stdout,*) 366 367 else 368 ! 369 ! geometry is set to the cell parameters read from stdin 370 ! 371 WRITE(stdout, '(3X,"ibrav = ",i4," cell parameters read from input file")') ibrav 372 373 h = at * alat 374 hold = h 375 376 end if 377 ! 378 ! generate true g-space 379 ! 380 call newinit( ht0%hmat, iverbosity = 1 ) 381 ! 382 CALL invmat( 3, h, ainv, deth ) 383 ! 384 344 format(3X,'ibrav = ',i4,' cell parameters ',/) 385 345 format(3(4x,f10.5)) 386 return 387 end subroutine init_geometry 388 389!----------------------------------------------------------------------- 390 391 subroutine newinit_x( h, iverbosity ) 392 ! 393 ! re-initialization of lattice parameters and g-space vectors. 394 ! Note that direct and reciprocal lattice primitive vectors 395 ! at, ainv, and corresponding quantities for small boxes 396 ! are recalculated according to the value of cell parameter h 397 ! 398 USE kinds, ONLY : DP 399 USE constants, ONLY : tpi 400 USE cell_base, ONLY : at, bg, omega, alat, tpiba2, & 401 cell_base_reinit 402 USE gvecw, ONLY : g2kin_init 403 USE gvect, ONLY : g, gg, ngm, mill 404 USE fft_base, ONLY : dfftp, dfftb 405 USE small_box, ONLY : small_box_set 406 USE smallbox_subs, ONLY : gcalb 407 USE io_global, ONLY : stdout, ionode 408 ! 409 implicit none 410 ! 411 REAL(DP), INTENT(IN) :: h(3,3) 412 INTEGER, INTENT(IN) :: iverbosity 413 ! 414 REAL(DP) :: rat1, rat2, rat3 415 INTEGER :: ig 416 ! 417 !WRITE( stdout, "(4x,'h from newinit')" ) 418 !do i=1,3 419 ! WRITE( stdout, '(3(4x,f12.7)' ) (h(i,j),j=1,3) 420 !enddo 421 ! 422 ! re-initialize the cell base module with the new geometry 423 ! 424 CALL cell_base_reinit( TRANSPOSE( h ) ) 425 ! 426 ! re-calculate G-vectors and kinetic energy 427 ! 428 do ig = 1, dfftp%ngm 429 g(:,ig)= mill(1,ig)*bg(:,1) + mill(2,ig)*bg(:,2) + mill(3,ig)*bg(:,3) 430 gg(ig)=g(1,ig)**2 + g(2,ig)**2 + g(3,ig)**2 431 enddo 432 ! 433 call g2kin_init ( gg, tpiba2 ) 434 ! 435 IF ( dfftb%nr1 == 0 .OR. dfftb%nr2 == 0 .OR. dfftb%nr3 == 0 ) RETURN 436 ! 437 ! generation of little box g-vectors 438 ! 439 rat1 = DBLE( dfftb%nr1 ) / DBLE( dfftp%nr1 ) 440 rat2 = DBLE( dfftb%nr2 ) / DBLE( dfftp%nr2 ) 441 rat3 = DBLE( dfftb%nr3 ) / DBLE( dfftp%nr3 ) 442 CALL small_box_set( alat, omega, at, rat1, rat2, rat3, tprint = ( iverbosity > 0 ) ) 443 ! 444 call gcalb ( ) 445 ! 446 ! pass new cell parameters to plugins 447 ! 448 CALL plugin_init_cell( ) 449 ! 450 return 451 end subroutine newinit_x 452 453 SUBROUTINE realspace_grids_info ( dfftp, dffts ) 454 455 ! Print info on local and global dimensions for real space grids 456 457 USE fft_types, ONLY: fft_type_descriptor 458 use io_global, only: stdout, ionode 459 USE fft_helper_subroutines, ONLY: fft_dist_info 460 461 IMPLICIT NONE 462 463 TYPE(fft_type_descriptor), INTENT(IN) :: dfftp, dffts 464 465 IF(ionode) THEN 466 WRITE( stdout,*) 467 WRITE( stdout,*) ' Real Mesh' 468 WRITE( stdout,*) ' ---------' 469 CALL fft_dist_info( dfftp, stdout ) 470 WRITE( stdout,*) 471 WRITE( stdout,*) ' Smooth Real Mesh' 472 WRITE( stdout,*) ' ----------------' 473 CALL fft_dist_info( dffts, stdout ) 474 END IF 475 476 RETURN 477 END SUBROUTINE realspace_grids_info 478 479 480