1! 2! Copyright (C) 2002-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 MODULE cell_base 10!------------------------------------------------------------------------------! 11 12 USE kinds, ONLY : DP 13 USE constants, ONLY : pi, bohr_radius_angs 14 USE io_global, ONLY : stdout 15! 16 IMPLICIT NONE 17 SAVE 18 ! 19 ! ibrav: index of the bravais lattice (see latgen.f90) 20 INTEGER :: ibrav 21 ! celldm: old-style parameters of the simulation cell (se latgen.f90) 22 REAL(DP) :: celldm(6) = (/ 0.0_DP,0.0_DP,0.0_DP,0.0_DP,0.0_DP,0.0_DP /) 23 ! traditional crystallographic cell parameters (alpha=cosbc and so on) 24 25 REAL(DP) :: a, b, c, cosab, cosac, cosbc 26 ! format of input cell parameters: 27 ! 'alat','bohr','angstrom' 28 CHARACTER(len=80) :: cell_units 29 ! alat: lattice parameter - often used to scale quantities, or 30 ! in combination to other parameters/constants to define new units 31 REAL(DP) :: alat = 0.0_DP 32 ! omega: volume of the simulation cell 33 REAl(DP) :: omega = 0.0_DP 34 ! tpiba: 2 PI/alat, tpiba2=tpiba^2 35 REAL(DP) :: tpiba = 0.0_DP, tpiba2 = 0.0_DP 36 ! direct and reciprocal lattice primitive vectors 37 ! at(:,i) are the lattice vectors of the simulation cell, a_i, 38 ! in alat units: a_i(:) = at(:,i)/alat 39 ! bg(:,i) are the reciprocal lattice vectors, b_i, 40 ! in tpiba=2pi/alat units: b_i(:) = bg(:,i)/tpiba 41 REAL(DP) :: at(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) ) 42 REAL(DP) :: bg(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) ) 43 ! 44 ! parameters for reference cell 45 REAL(DP) :: ref_tpiba2 = 0.0_DP 46 REAL(DP) :: ref_at(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) ) 47 REAL(DP) :: ref_bg(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) ) 48 ! 49 ! parameter to store tpiba2 calculated from the input cell parameter 50 ! used in emass_preconditioning, required for restarting variable cell calculation correctly in CP 51 REAL(DP) :: init_tpiba2 = 0.0_DP 52 ! 53! ------------------------------------------------------------------------- 54! ... periodicity box 55! ... In the matrix "a" every row is the vector of each side of 56! ... the cell in the real space 57 58 TYPE boxdimensions 59 REAL(DP) :: a(3,3) ! direct lattice generators 60 REAL(DP) :: m1(3,3) ! reciprocal lattice generators 61 REAL(DP) :: omega ! cell volume = determinant of a 62 REAL(DP) :: g(3,3) ! metric tensor 63 REAL(DP) :: gvel(3,3) ! metric velocity 64 REAL(DP) :: pail(3,3) ! stress tensor ( scaled coor. ) 65 REAL(DP) :: paiu(3,3) ! stress tensor ( cartesian coor. ) 66 REAL(DP) :: hmat(3,3) ! cell parameters ( transpose of "a" ) 67 REAL(DP) :: hvel(3,3) ! cell velocity 68 REAL(DP) :: hinv(3,3) 69 REAL(DP) :: deth 70 INTEGER :: perd(3) 71 END TYPE boxdimensions 72 73 ! The following relations should always be kept valid: 74 ! h = at*alat; ainv = h^(-1); ht=transpose(h) 75 REAL(DP) :: h(3,3) = 0.0_DP ! simulation cell at time t 76 REAL(DP) :: ainv(3,3) = 0.0_DP 77 REAL(DP) :: hold(3,3) = 0.0_DP ! simulation cell at time t-delt 78 REAL(DP) :: hnew(3,3) = 0.0_DP ! simulation cell at time t+delt 79 REAL(DP) :: velh(3,3) = 0.0_DP ! simulation cell velocity 80 REAL(DP) :: deth = 0.0_DP ! determinant of h ( cell volume ) 81 82 INTEGER :: iforceh(3,3) = 1 ! if iforceh( i, j ) = 0 then h( i, j ) 83 ! is not allowed to move 84 LOGICAL :: enforce_ibrav = .FALSE.! True if ibrav representation is fix 85 LOGICAL :: fix_volume = .FALSE.! True if cell volume is kept fixed 86 LOGICAL :: fix_area = .FALSE. ! True if area in xy plane is kept constant 87 LOGICAL :: isotropic = .FALSE. ! True if volume option is chosen for cell_dofree 88 REAL(DP) :: wmass = 0.0_DP ! cell fictitious mass 89 REAL(DP) :: press = 0.0_DP ! external pressure 90 91 REAL(DP) :: frich = 0.0_DP ! friction parameter for cell damped dynamics 92 REAL(DP) :: greash = 1.0_DP ! greas parameter for damped dynamics 93 94 LOGICAL :: tcell_base_init = .FALSE. 95 96 INTERFACE cell_init 97 MODULE PROCEDURE cell_init_ht, cell_init_a 98 END INTERFACE 99 100 INTERFACE pbcs 101 MODULE PROCEDURE pbcs_components, pbcs_vectors 102 END INTERFACE 103 104 INTERFACE s_to_r 105 MODULE PROCEDURE s_to_r1, s_to_r1b, s_to_r3 106 END INTERFACE 107 108 INTERFACE r_to_s 109 MODULE PROCEDURE r_to_s1, r_to_s1b, r_to_s3 110 END INTERFACE 111!------------------------------------------------------------------------------! 112 CONTAINS 113!------------------------------------------------------------------------------! 114! 115 SUBROUTINE cell_base_init( ibrav_, celldm_, a_, b_, c_, cosab_, cosac_, & 116 cosbc_, trd_ht, rd_ht, cell_units_ ) 117 ! 118 ! ... initialize cell_base module variables, set up crystal lattice 119 ! 120 121 IMPLICIT NONE 122 INTEGER, INTENT(IN) :: ibrav_ 123 REAL(DP), INTENT(IN) :: celldm_ (6) 124 LOGICAL, INTENT(IN) :: trd_ht 125 REAL(DP), INTENT(IN) :: rd_ht (3,3) 126 CHARACTER(LEN=*), INTENT(IN) :: cell_units_ 127 REAL(DP), INTENT(IN) :: a_ , b_ , c_ , cosab_, cosac_, cosbc_ 128 129 REAL(DP) :: units 130 ! 131 IF ( ibrav_ == 0 .and. .not. trd_ht ) THEN 132 CALL errore('cell_base_init', 'ibrav=0: must read cell parameters', 1) 133 ELSE IF ( ibrav_ /= 0 .and. trd_ht ) THEN 134 CALL errore('cell_base_init', 'redundant data for cell parameters', 2) 135 END IF 136 ! 137 ibrav = ibrav_ 138 celldm = celldm_ 139 a = a_ ; b = b_ ; c = c_ ; cosab = cosab_ ; cosac = cosac_ ; cosbc = cosbc_ 140 cell_units = cell_units_ 141 units = 0.0_DP 142 ! 143 IF ( trd_ht ) THEN 144 ! 145 ! ... crystal lattice vectors read from input: find units 146 ! 147 SELECT CASE ( TRIM( cell_units ) ) 148 CASE ( 'bohr' ) 149 IF( celldm( 1 ) /= 0.0_DP .OR. a /= 0.0_dp ) CALL errore & 150 ('cell_base_init','lattice parameter specified twice',1) 151 units = 1.0_DP 152 CASE ( 'angstrom' ) 153 IF( celldm( 1 ) /= 0.0_DP .OR. a /= 0.0_dp ) CALL errore & 154 ('cell_base_init','lattice parameter specified twice',2) 155 units = 1.0_DP / bohr_radius_angs 156 CASE ( 'alat' ) 157 IF( celldm( 1 ) /= 0.0_DP ) THEN 158 units = celldm( 1 ) 159 ELSE IF ( a /= 0.0_dp ) THEN 160 units = a / bohr_radius_angs 161 ELSE 162 CALL errore ('cell_base_init', & 163 'lattice parameter not specified',1) 164 END IF 165 ! following case is deprecated and should be removed 166 CASE ( 'none' ) 167 ! cell_units is 'none' if nothing was specified 168 IF( celldm( 1 ) /= 0.0_DP ) THEN 169 units = celldm( 1 ) 170 cell_units = 'alat' 171 ELSE IF ( a /= 0.0_dp ) THEN 172 units = a / bohr_radius_angs 173 cell_units = 'alat' 174 ELSE 175 units = 1.0_DP 176 cell_units = 'bohr' 177 END IF 178 ! 179 CASE DEFAULT 180 CALL errore ('cell_base_init', & 181 'unexpected cell_units '//TRIM(cell_units),1) 182 END SELECT 183 ! 184 ! ... Beware the transpose operation between matrices ht and at! 185 ! 186 at = TRANSPOSE( rd_ht ) * units 187 ! 188 ! ... at is in atomic units: find alat, bring at to alat units, find omega 189 ! 190 IF( celldm( 1 ) /= 0.0_DP ) THEN 191 alat = celldm( 1 ) 192 ELSE IF ( a /= 0.0_dp ) THEN 193 alat = a / bohr_radius_angs 194 ELSE 195 alat = SQRT ( at(1,1)**2+at(2,1)**2+at(3,1)**2 ) 196 END IF 197 ! for compatibility: celldm still used in phonon etc 198 celldm(1) = alat 199 ! 200 at(:,:) = at(:,:) / alat 201 CALL volume( alat, at(1,1), at(1,2), at(1,3), omega ) 202 ! 203 ELSE 204 ! 205 ! ... crystal lattice vectors via ibrav + celldm parameters 206 ! 207 IF ( celldm(1) == 0.D0 .and. a /= 0.D0 ) THEN 208 ! 209 ! ... convert crystallographic parameters into celldm parameters 210 ! 211 CALL abc2celldm ( ibrav, a,b,c,cosab,cosac,cosbc, celldm ) 212 ! 213 ELSE IF ( celldm(1) /= 0.D0 .and. a /= 0.D0 ) THEN 214 ! 215 CALL errore( 'input', 'do not specify both celldm and a,b,c!', 1 ) 216 ! 217 END IF 218 ! 219 ! ... generate at (in atomic units) from ibrav and celldm 220 ! 221 CALL latgen( ibrav, celldm, at(1,1), at(1,2), at(1,3), omega ) 222 ! 223 ! ... define lattice constants alat, divide at by alat 224 ! 225 alat = celldm(1) 226 at(:,:) = at(:,:) / alat 227 ! 228 END IF 229 IF ( alat < 1.9_dp ) CALL infomsg ('cell_base_init', & 230 'DEPRECATED: use true lattice parameter, not A to a.u. conversion factor') 231 ! 232 ! ... Generate the reciprocal lattice vectors 233 ! 234 CALL recips( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) ) 235 ! 236 tpiba = 2.0_DP * pi / alat 237 tpiba2 = tpiba * tpiba 238 init_tpiba2 = tpiba2 ! BS : this is used in CPV/src/init_run.f90 239 RETURN 240 ! 241 END SUBROUTINE cell_base_init 242 ! 243 SUBROUTINE ref_cell_base_init( ref_alat, rd_ref_ht, ref_cell_units ) 244 ! 245 ! ... initialize cell_base module variables, set up crystal lattice 246 ! 247 248 IMPLICIT NONE 249 REAL(DP), INTENT(IN) :: rd_ref_ht (3,3) 250 REAL(DP), INTENT(INOUT) :: ref_alat 251 CHARACTER(LEN=*), INTENT(IN) :: ref_cell_units 252 253 REAL(DP) :: units, ref_omega 254 ! 255 ! ... reference cell lattice vectors read from REF_CELL_PARAMETERS Card: find units 256 ! 257 SELECT CASE ( TRIM( ref_cell_units ) ) 258 ! 259 CASE ( 'bohr' ) 260 units = 1.0_DP 261 CASE ( 'angstrom' ) 262 units = 1.0_DP / bohr_radius_angs 263 CASE DEFAULT 264 IF( ref_alat .GT. 0.0_DP ) THEN 265 units = ref_alat 266 ELSE 267 CALL errore('ref_cell_base_init', 'ref_alat must be set to a positive value (in A.U.) in SYSTEM namelist', 1) 268 END IF 269 ! 270 END SELECT 271 ! 272 ! ... Beware the transpose operation between matrices ht and at! 273 ! 274 ref_at = TRANSPOSE( rd_ref_ht ) * units 275 ! 276 ! ... ref_at is in atomic units: find ref_alat, bring ref_at to ref_alat units 277 ! 278 ref_alat = SQRT ( ref_at(1,1)**2+ref_at(2,1)**2+ref_at(3,1)**2 ) 279 ! 280 ref_at(:,:) = ref_at(:,:) / ref_alat 281 ! 282 ! ... Generate the reciprocal lattice vectors from the reference cell 283 ! 284 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) ) 285 ! 286 ref_tpiba2 = (2.0_DP * pi / ref_alat)**2 287 ! 288 CALL volume( ref_alat, ref_at(1,1), ref_at(1,2), ref_at(1,3), ref_omega ) 289 ! 290 WRITE( stdout, * ) 291 WRITE( stdout, '(3X,"Reference Cell read from REF_CELL_PARAMETERS Card")' ) 292 WRITE( stdout, '(3X,"Reference Cell alat =",F14.8,1X,"A.U.")' ) ref_alat 293 WRITE( stdout, '(3X,"ref_cell_a1 = ",1X,3f14.8)' ) ref_at(:,1)*ref_alat 294 WRITE( stdout, '(3X,"ref_cell_a2 = ",1X,3f14.8)' ) ref_at(:,2)*ref_alat 295 WRITE( stdout, '(3X,"ref_cell_a3 = ",1X,3f14.8)' ) ref_at(:,3)*ref_alat 296 WRITE( stdout, * ) 297 WRITE( stdout, '(3X,"ref_cell_b1 = ",1X,3f14.8)' ) ref_bg(:,1)/ref_alat 298 WRITE( stdout, '(3X,"ref_cell_b2 = ",1X,3f14.8)' ) ref_bg(:,2)/ref_alat 299 WRITE( stdout, '(3X,"ref_cell_b3 = ",1X,3f14.8)' ) ref_bg(:,3)/ref_alat 300 WRITE( stdout, '(3X,"Reference Cell Volume",F16.8,1X,"A.U.")' ) ref_omega 301 ! 302 RETURN 303 ! 304 END SUBROUTINE ref_cell_base_init 305!------------------------------------------------------------------------------! 306! ... set box 307! ... box%m1(i,1) == b1(i) COLUMN are B vectors 308! ... box%a(1,i) == a1(i) ROW are A vector 309! ... box%omega == volume 310! ... box%g(i,j) == metric tensor G 311!------------------------------------------------------------------------------! 312 313 SUBROUTINE cell_init_ht( what, box, hval ) 314 TYPE (boxdimensions) :: box 315 REAL(DP), INTENT(IN) :: hval(3,3) 316 CHARACTER, INTENT(IN) :: what 317 IF( what == 't' .OR. what == 'T' ) THEN 318 ! hval == ht 319 box%a = hval 320 box%hmat = TRANSPOSE( hval ) 321 ELSE 322 ! hval == hmat 323 box%hmat = hval 324 box%a = TRANSPOSE( hval ) 325 END IF 326 CALL gethinv( box ) 327 box%g = MATMUL( box%a(:,:), box%hmat(:,:) ) 328 box%gvel = 0.0_DP 329 box%hvel = 0.0_DP 330 box%pail = 0.0_DP 331 box%paiu = 0.0_DP 332 RETURN 333 END SUBROUTINE cell_init_ht 334 335!------------------------------------------------------------------------------! 336 337 SUBROUTINE cell_init_a( alat, at, box ) 338 TYPE (boxdimensions) :: box 339 REAL(DP), INTENT(IN) :: alat, at(3,3) 340 INTEGER :: i 341 DO i=1,3 342 ! this is HT: the rows are the lattice vectors 343 box%a(1,i) = at(i,1)*alat 344 box%a(2,i) = at(i,2)*alat 345 box%a(3,i) = at(i,3)*alat 346 ! this is H : the column are the lattice vectors 347 box%hmat(i,1) = at(i,1)*alat 348 box%hmat(i,2) = at(i,2)*alat 349 box%hmat(i,3) = at(i,3)*alat 350 END DO 351 box%pail = 0.0_DP 352 box%paiu = 0.0_DP 353 box%hvel = 0.0_DP 354 CALL gethinv(box) 355 box%g = MATMUL( box%a(:,:), box%hmat(:,:) ) 356 box%gvel = 0.0_DP 357 RETURN 358 END SUBROUTINE cell_init_a 359 360!------------------------------------------------------------------------------! 361 362 SUBROUTINE r_to_s1 (r,s,box) 363 REAL(DP), intent(out) :: S(3) 364 REAL(DP), intent(in) :: R(3) 365 type (boxdimensions), intent(in) :: box 366 integer i,j 367 DO I=1,3 368 S(I) = 0.0_DP 369 DO J=1,3 370 S(I) = S(I) + R(J)*box%m1(J,I) 371 END DO 372 END DO 373 RETURN 374 END SUBROUTINE r_to_s1 375 376!------------------------------------------------------------------------------! 377 378 SUBROUTINE r_to_s3 ( r, s, nat, hinv ) 379 REAL(DP), intent(out) :: S(:,:) 380 INTEGER, intent(in) :: nat 381 REAL(DP), intent(in) :: R(:,:) 382 REAL(DP), intent(in) :: hinv(:,:) ! hinv = TRANSPOSE( box%m1 ) 383 integer :: i, j, ia 384 DO ia = 1, nat 385 DO i=1,3 386 S(i,ia) = 0.0_DP 387 DO j=1,3 388 S(i,ia) = S(i,ia) + R(j,ia)*hinv(i,j) 389 END DO 390 END DO 391 END DO 392 RETURN 393 END SUBROUTINE r_to_s3 394 395!------------------------------------------------------------------------------! 396 397 SUBROUTINE r_to_s1b ( r, s, hinv ) 398 REAL(DP), intent(out) :: S(:) 399 REAL(DP), intent(in) :: R(:) 400 REAL(DP), intent(in) :: hinv(:,:) ! hinv = TRANSPOSE( box%m1 ) 401 integer :: i, j 402 DO I=1,3 403 S(I) = 0.0_DP 404 DO J=1,3 405 S(I) = S(I) + R(J)*hinv(i,j) 406 END DO 407 END DO 408 RETURN 409 END SUBROUTINE r_to_s1b 410 411 412!------------------------------------------------------------------------------! 413 414 SUBROUTINE s_to_r1 (S,R,box) 415 REAL(DP), intent(in) :: S(3) 416 REAL(DP), intent(out) :: R(3) 417 type (boxdimensions), intent(in) :: box 418 integer i,j 419 DO I=1,3 420 R(I) = 0.0_DP 421 DO J=1,3 422 R(I) = R(I) + S(J)*box%a(J,I) 423 END DO 424 END DO 425 RETURN 426 END SUBROUTINE s_to_r1 427 428!------------------------------------------------------------------------------! 429 430 SUBROUTINE s_to_r1b (S,R,h) 431 REAL(DP), intent(in) :: S(3) 432 REAL(DP), intent(out) :: R(3) 433 REAL(DP), intent(in) :: h(:,:) ! h = TRANSPOSE( box%a ) 434 integer i,j 435 DO I=1,3 436 R(I) = 0.0_DP 437 DO J=1,3 438 R(I) = R(I) + S(J)*h(I,j) 439 END DO 440 END DO 441 RETURN 442 END SUBROUTINE s_to_r1b 443 444!------------------------------------------------------------------------------! 445 446 SUBROUTINE s_to_r3 ( S, R, nat, h ) 447 REAL(DP), intent(in) :: S(:,:) 448 INTEGER, intent(in) :: nat 449 REAL(DP), intent(out) :: R(:,:) 450 REAL(DP), intent(in) :: h(:,:) ! h = TRANSPOSE( box%a ) 451 integer :: i, j, ia 452 DO ia = 1, nat 453 DO I = 1, 3 454 R(I,ia) = 0.0_DP 455 DO J = 1, 3 456 R(I,ia) = R(I,ia) + S(J,ia) * h(I,j) 457 END DO 458 END DO 459 END DO 460 RETURN 461 END SUBROUTINE s_to_r3 462 463 464! 465!------------------------------------------------------------------------------! 466! 467 468 SUBROUTINE gethinv(box) 469 USE matrix_inversion 470 IMPLICIT NONE 471 TYPE (boxdimensions), INTENT (INOUT) :: box 472 ! 473 CALL invmat( 3, box%a, box%m1, box%omega ) 474 box%deth = box%omega 475 box%hinv = TRANSPOSE( box%m1 ) 476 ! 477 RETURN 478 END SUBROUTINE gethinv 479 480 481 FUNCTION get_volume( hmat ) 482 IMPLICIT NONE 483 REAL(DP) :: get_volume 484 REAL(DP) :: hmat( 3, 3 ) 485 get_volume = hmat(1,1)*(hmat(2,2)*hmat(3,3)-hmat(2,3)*hmat(3,2)) + & 486 hmat(1,2)*(hmat(2,3)*hmat(3,1)-hmat(2,1)*hmat(3,3)) + & 487 hmat(1,3)*(hmat(2,1)*hmat(3,2)-hmat(2,2)*hmat(3,1)) 488 RETURN 489 END FUNCTION get_volume 490! 491!------------------------------------------------------------------------------! 492! 493 FUNCTION pbc(rin,box,nl) RESULT (rout) 494 IMPLICIT NONE 495 TYPE (boxdimensions) :: box 496 REAL (DP) :: rin(3) 497 REAL (DP) :: rout(3), s(3) 498 INTEGER, OPTIONAL :: nl(3) 499 500 s = matmul(box%hinv(:,:),rin) 501 s = s - box%perd*nint(s) 502 rout = matmul(box%hmat(:,:),s) 503 IF (present(nl)) THEN 504 s = REAL( nl, DP ) 505 rout = rout + matmul(box%hmat(:,:),s) 506 END IF 507 END FUNCTION pbc 508! 509!------------------------------------------------------------------------------! 510! 511 SUBROUTINE get_cell_param(box,cell,ang) 512 IMPLICIT NONE 513 TYPE(boxdimensions), INTENT(in) :: box 514 REAL(DP), INTENT(out), DIMENSION(3) :: cell 515 REAL(DP), INTENT(out), DIMENSION(3), OPTIONAL :: ang 516! This code gets the cell parameters given the h-matrix: 517! a 518 cell(1)=sqrt(box%hmat(1,1)*box%hmat(1,1)+box%hmat(2,1)*box%hmat(2,1) & 519 +box%hmat(3,1)*box%hmat(3,1)) 520! b 521 cell(2)=sqrt(box%hmat(1,2)*box%hmat(1,2)+box%hmat(2,2)*box%hmat(2,2) & 522 +box%hmat(3,2)*box%hmat(3,2)) 523! c 524 cell(3)=sqrt(box%hmat(1,3)*box%hmat(1,3)+box%hmat(2,3)*box%hmat(2,3) & 525 +box%hmat(3,3)*box%hmat(3,3)) 526 IF (PRESENT(ang)) THEN 527! gamma 528 ang(1)=acos((box%hmat(1,1)*box%hmat(1,2)+ & 529 box%hmat(2,1)*box%hmat(2,2) & 530 +box%hmat(3,1)*box%hmat(3,2))/(cell(1)*cell(2))) 531! beta 532 ang(2)=acos((box%hmat(1,1)*box%hmat(1,3)+ & 533 box%hmat(2,1)*box%hmat(2,3) & 534 +box%hmat(3,1)*box%hmat(3,3))/(cell(1)*cell(3))) 535! alpha 536 ang(3)=acos((box%hmat(1,2)*box%hmat(1,3)+ & 537 box%hmat(2,2)*box%hmat(2,3) & 538 +box%hmat(3,2)*box%hmat(3,3))/(cell(2)*cell(3))) 539! ang=ang*180.0_DP/pi 540 541 ENDIF 542 END SUBROUTINE get_cell_param 543 544!------------------------------------------------------------------------------! 545 546 SUBROUTINE pbcs_components(x1, y1, z1, x2, y2, z2, m) 547! ... This subroutine compute the periodic boundary conditions in the scaled 548! ... variables system 549 USE kinds 550 INTEGER, INTENT(IN) :: M 551 REAL(DP), INTENT(IN) :: X1,Y1,Z1 552 REAL(DP), INTENT(OUT) :: X2,Y2,Z2 553 REAL(DP) MIC 554 MIC = REAL( M, DP ) 555 X2 = X1 - DNINT(X1/MIC)*MIC 556 Y2 = Y1 - DNINT(Y1/MIC)*MIC 557 Z2 = Z1 - DNINT(Z1/MIC)*MIC 558 RETURN 559 END SUBROUTINE pbcs_components 560 561!------------------------------------------------------------------------------! 562 563 SUBROUTINE pbcs_vectors(v, w, m) 564! ... This subroutine compute the periodic boundary conditions in the scaled 565! ... variables system 566 USE kinds 567 INTEGER, INTENT(IN) :: m 568 REAL(DP), INTENT(IN) :: v(3) 569 REAL(DP), INTENT(OUT) :: w(3) 570 REAL(DP) :: MIC 571 MIC = REAL( M, DP ) 572 w(1) = v(1) - DNINT(v(1)/MIC)*MIC 573 w(2) = v(2) - DNINT(v(2)/MIC)*MIC 574 w(3) = v(3) - DNINT(v(3)/MIC)*MIC 575 RETURN 576 END SUBROUTINE pbcs_vectors 577 578!------------------------------------------------------------------------------! 579 580 SUBROUTINE set_h_ainv() 581 ! 582 ! CP-PW compatibility: align CP arrays H and ainv to at and bg 583 ! 584 IMPLICIT NONE 585 ! 586 !write(stdout,*) 'alat=',alat 587 !write(stdout,*) 'at=',at 588 !write(stdout,*) 'bg=',bg 589 ! 590 h(:,:) = at(:,:)*alat 591 ! 592 ainv(1,:) = bg(:,1)/alat 593 ainv(2,:) = bg(:,2)/alat 594 ainv(3,:) = bg(:,3)/alat 595 ! 596 END SUBROUTINE set_h_ainv 597 598!------------------------------------------------------------------------------! 599 600 SUBROUTINE cell_dyn_init( trd_ht, rd_ht, wc_ , total_ions_mass , press_ , & 601 frich_ , greash_ , cell_dofree ) 602 603 USE constants, ONLY: au_gpa, amu_au 604 USE io_global, ONLY: stdout 605 606 IMPLICIT NONE 607 CHARACTER(LEN=*), INTENT(IN) :: cell_dofree 608 LOGICAL, INTENT(IN) :: trd_ht 609 REAL(DP), INTENT(IN) :: rd_ht (3,3) 610 REAL(DP), INTENT(IN) :: wc_ , frich_ , greash_ , total_ions_mass 611 REAL(DP), INTENT(IN) :: press_ ! external pressure from input 612 ! ( in KBar = 0.1 GPa ) 613 INTEGER :: j 614 ! 615 press = press_ / 10.0_DP ! convert press in KBar to GPa 616 press = press / au_gpa ! convert to AU 617 ! frich = frich_ ! for the time being this is set elsewhere 618 greash = greash_ 619 620 WRITE( stdout, 105 ) 621 WRITE( stdout, 110 ) press_ 622105 format(/,3X,'Simulation Cell Parameters (from input)') 623110 format( 3X,'external pressure = ',f15.2,' [KBar]') 624 625 wmass = wc_ 626 IF( wmass == 0.0_DP ) THEN 627 wmass = 3.0_DP / (4.0_DP * pi**2 ) * total_ions_mass 628 wmass = wmass * AMU_AU 629 WRITE( stdout,130) wmass 630 ELSE 631 WRITE( stdout,120) wmass 632 END IF 633120 format(3X,'wmass (read from input) = ',f15.2,' [AU]') 634130 format(3X,'wmass (calculated) = ',f15.2,' [AU]') 635 636 IF( wmass <= 0.0_DP ) & 637 CALL errore(' cell_dyn_init',' wmass out of range ',0) 638 639 IF ( trd_ht ) THEN 640 ! 641 WRITE( stdout, 210 ) 642 WRITE( stdout, 220 ) ( rd_ht( 1, j ), j = 1, 3 ) 643 WRITE( stdout, 220 ) ( rd_ht( 2, j ), j = 1, 3 ) 644 WRITE( stdout, 220 ) ( rd_ht( 3, j ), j = 1, 3 ) 645 ! 646210 format(3X,'initial cell from CELL_PARAMETERS card') 647220 format(3X,3F14.8) 648 ! 649 END IF 650 ! 651 ainv(1,:) = bg(:,1)/alat 652 ainv(2,:) = bg(:,2)/alat 653 ainv(3,:) = bg(:,3)/alat 654 ! 655 CALL init_dofree ( cell_dofree ) 656 ! 657 tcell_base_init = .TRUE. 658 659 WRITE( stdout, 300 ) ibrav 660 WRITE( stdout, 305 ) alat 661 WRITE( stdout, 310 ) at(:,1)*alat 662 WRITE( stdout, 320 ) at(:,2)*alat 663 WRITE( stdout, 330 ) at(:,3)*alat 664 WRITE( stdout, * ) 665 WRITE( stdout, 350 ) bg(:,1)/alat 666 WRITE( stdout, 360 ) bg(:,2)/alat 667 WRITE( stdout, 370 ) bg(:,3)/alat 668 WRITE( stdout, 340 ) omega 669300 FORMAT( 3X, 'ibrav = ',I4) 670305 FORMAT( 3X, 'alat = ',F14.8) 671310 FORMAT( 3X, 'a1 = ',3F14.8) 672320 FORMAT( 3X, 'a2 = ',3F14.8) 673330 FORMAT( 3X, 'a3 = ',3F14.8) 674350 FORMAT( 3X, 'b1 = ',3F14.8) 675360 FORMAT( 3X, 'b2 = ',3F14.8) 676370 FORMAT( 3X, 'b3 = ',3F14.8) 677340 FORMAT( 3X, 'omega = ',F16.8) 678 679 680 RETURN 681 END SUBROUTINE cell_dyn_init 682 683!------------------------------------------------------------------------------! 684 SUBROUTINE init_dofree ( cell_dofree ) 685 686 ! set constraints on cell dynamics/optimization 687 688 CHARACTER(LEN=*), INTENT(IN) :: cell_dofree 689 690 SELECT CASE ( TRIM( cell_dofree ) ) 691 692 CASE ( 'all', 'default' ) 693 iforceh = 1 694 CASE ( 'ibrav') 695 iforceh = 1 696 enforce_ibrav = .true. 697 CASE ( 'shape' ) 698 iforceh = 1 699 fix_volume = .true. 700! 2DSHAPE: CASE FOR SHAPE CHANGE IN xy PLANE WITH CONST AREA 701! contribution from Richard Charles Andrew 702! Physics Department, University of Pretoria 703! South Africa, august 2012. 704 CASE ( '2Dshape' ) 705 iforceh = 1 706 iforceh(3,3) = 0 707 iforceh(1,3) = 0 708 iforceh(3,1) = 0 709 iforceh(2,3) = 0 710 iforceh(3,2) = 0 711 fix_area = .true. 712! 2DSHAPE 713 CASE ( 'volume' ) 714 !CALL errore(' init_dofree ', & 715 ! ' cell_dofree = '//TRIM(cell_dofree)//' not yet implemented ', 1 ) 716 IF ( ibrav /= 1 ) THEN 717 CALL errore('cell_dofree', 'Isotropic expansion is only allowed for ibrav=1; i.e. for simple cubic', 1) 718 END IF 719 iforceh = 0 720 iforceh(1,1) = 1 721 iforceh(2,2) = 1 722 iforceh(3,3) = 1 723 isotropic = .TRUE. 724 CASE ('x') 725 iforceh = 0 726 iforceh(1,1) = 1 727 CASE ('y') 728 iforceh = 0 729 iforceh(2,2) = 1 730 CASE ('z') 731 iforceh = 0 732 iforceh(3,3) = 1 733 CASE ('xy') 734 iforceh = 0 735 iforceh(1,1) = 1 736 iforceh(2,2) = 1 737 ! ... if you want the entire xy plane to be free, uncomment: 738 ! iforceh(1,2) = 1 739 ! iforceh(2,1) = 1 740! 2DSHAPE THE ENTIRE xy PLANE IS FREE 741 CASE ('2Dxy') 742 iforceh = 0 743 iforceh(1,1) = 1 744 iforceh(2,2) = 1 745 iforceh(1,2) = 1 746 iforceh(2,1) = 1 747! 2DSHAPE 748 CASE ('xz') 749 iforceh = 0 750 iforceh(1,1) = 1 751 iforceh(3,3) = 1 752 CASE ('yz') 753 iforceh = 0 754 iforceh(2,2) = 1 755 iforceh(3,3) = 1 756 CASE ('xyz') 757 iforceh = 0 758 iforceh(1,1) = 1 759 iforceh(2,2) = 1 760 iforceh(3,3) = 1 761! epitaxial constraints (2 axes fixed, one free) 762! added by ulrich.aschauer@dcb.unibe.ch on 2018-02-02 763 CASE ('epitaxial_ab') 764 !fix the a and b axis while allowing c to change 765 iforceh = 0 766 iforceh(1,3) = 1 767 iforceh(2,3) = 1 768 iforceh(3,3) = 1 769 CASE ('epitaxial_ac') 770 !fix the a and c axis while allowing b to change 771 iforceh = 0 772 iforceh(1,2) = 1 773 iforceh(2,2) = 1 774 iforceh(3,2) = 1 775 CASE ('epitaxial_bc') 776 !fix the b and c axis while allowing a to change 777 iforceh = 0 778 iforceh(1,1) = 1 779 iforceh(2,1) = 1 780 iforceh(3,1) = 1 781 CASE DEFAULT 782 CALL errore(' init_dofree ',' unknown cell_dofree '//TRIM(cell_dofree), 1 ) 783 784 END SELECT 785 END SUBROUTINE init_dofree 786 787!------------------------------------------------------------------------------! 788 789 SUBROUTINE cell_base_reinit( ht ) 790 791 USE control_flags, ONLY: iverbosity 792 793 IMPLICIT NONE 794 REAL(DP), INTENT(IN) :: ht (3,3) 795 796 INTEGER :: j 797 798 alat = sqrt( ht(1,1)*ht(1,1) + ht(1,2)*ht(1,2) + ht(1,3)*ht(1,3) ) 799 tpiba = 2.0_DP * pi / alat 800 tpiba2 = tpiba * tpiba 801 ! 802 IF( iverbosity > 2 ) THEN 803 WRITE( stdout, 210 ) 804 WRITE( stdout, 220 ) ( ht( 1, j ), j = 1, 3 ) 805 WRITE( stdout, 220 ) ( ht( 2, j ), j = 1, 3 ) 806 WRITE( stdout, 220 ) ( ht( 3, j ), j = 1, 3 ) 807 END IF 808210 format(3X,'Simulation cell parameters with the new cell:') 809220 format(3X,3F14.8) 810 811 ! matrix "ht" used in CP is the transpose of matrix "at" 812 ! times the lattice parameter "alat"; matrix "ainv" is "bg" divided alat 813 ! 814 at = TRANSPOSE( ht ) / alat 815 ! 816 CALL recips( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) ) 817 CALL volume( alat, at(1,1), at(1,2), at(1,3), deth ) 818 omega = deth 819 ! 820 ainv(1,:) = bg(:,1)/alat 821 ainv(2,:) = bg(:,2)/alat 822 ainv(3,:) = bg(:,3)/alat 823 ! 824 IF( iverbosity > 2 ) THEN 825 WRITE( stdout, 305 ) alat 826 WRITE( stdout, 310 ) at(:,1)*alat 827 WRITE( stdout, 320 ) at(:,2)*alat 828 WRITE( stdout, 330 ) at(:,3)*alat 829 WRITE( stdout, * ) 830 WRITE( stdout, 350 ) bg(:,1)/alat 831 WRITE( stdout, 360 ) bg(:,2)/alat 832 WRITE( stdout, 370 ) bg(:,3)/alat 833 WRITE( stdout, 340 ) omega 834 END IF 835 836305 FORMAT( 3X, 'alat = ',F14.8) 837310 FORMAT( 3X, 'a1 = ',3F14.8) 838320 FORMAT( 3X, 'a2 = ',3F14.8) 839330 FORMAT( 3X, 'a3 = ',3F14.8) 840350 FORMAT( 3X, 'b1 = ',3F14.8) 841360 FORMAT( 3X, 'b2 = ',3F14.8) 842370 FORMAT( 3X, 'b3 = ',3F14.8) 843340 FORMAT( 3X, 'omega = ',F14.8) 844 845 846 RETURN 847 END SUBROUTINE cell_base_reinit 848 849 850 851!------------------------------------------------------------------------------! 852 853 SUBROUTINE cell_steepest( hnew, h, delt, iforceh, fcell ) 854 REAL(DP), INTENT(OUT) :: hnew(3,3) 855 REAL(DP), INTENT(IN) :: h(3,3), fcell(3,3) 856 INTEGER, INTENT(IN) :: iforceh(3,3) 857 REAL(DP), INTENT(IN) :: delt 858 INTEGER :: i, j 859 REAL(DP) :: dt2,fiso 860 dt2 = delt * delt 861 ! 862 IF( isotropic ) THEN 863 ! 864 ! Isotropic force on the cell 865 ! 866 fiso = (fcell(1,1)+fcell(2,2)+fcell(3,3))/3.0_DP 867 ! 868 DO j=1,3 869 DO i=1,3 870 hnew(i,j) = h(i,j) + dt2 * fiso * REAL( iforceh(i,j), DP ) 871 ENDDO 872 ENDDO 873 ! 874 ELSE 875 ! 876 DO j=1,3 877 DO i=1,3 878 hnew(i,j) = h(i,j) + dt2 * fcell(i,j) * REAL( iforceh(i,j), DP ) 879 ENDDO 880 ENDDO 881 ! 882 END IF 883 ! 884 RETURN 885 END SUBROUTINE cell_steepest 886 887 888!------------------------------------------------------------------------------! 889 890 SUBROUTINE cell_verlet( hnew, h, hold, delt, iforceh, fcell, frich, tnoseh, hnos ) 891 REAL(DP), INTENT(OUT) :: hnew(3,3) 892 REAL(DP), INTENT(IN) :: h(3,3), hold(3,3), hnos(3,3), fcell(3,3) 893 INTEGER, INTENT(IN) :: iforceh(3,3) 894 REAL(DP), INTENT(IN) :: frich, delt 895 LOGICAL, INTENT(IN) :: tnoseh 896 897 REAL(DP) :: htmp(3,3) 898 REAL(DP) :: verl1, verl2, verl3, dt2, ftmp, v1, v2, v3, fiso 899 INTEGER :: i, j 900 901 dt2 = delt * delt 902 903 IF( tnoseh ) THEN 904 ftmp = 0.0_DP 905 htmp = hnos 906 ELSE 907 ftmp = frich 908 htmp = 0.0_DP 909 END IF 910 911 verl1 = 2.0_DP / ( 1.0_DP + ftmp ) 912 verl2 = 1.0_DP - verl1 913 verl3 = dt2 / ( 1.0_DP + ftmp ) 914 verl1 = verl1 - 1.0_DP 915 916 IF( isotropic ) THEN 917 ! 918 fiso = (fcell(1,1)+fcell(2,2)+fcell(3,3))/3.0_DP 919 ! 920 DO j=1,3 921 DO i=1,3 922 v1 = verl1 * h(i,j) 923 v2 = verl2 * hold(i,j) 924 v3 = verl3 * ( fiso - htmp(i,j) ) 925 hnew(i,j) = h(i,j) + ( v1 + v2 + v3 ) * REAL( iforceh(i,j), DP ) 926 ENDDO 927 ENDDO 928 ! 929 ELSE 930 ! 931 DO j=1,3 932 DO i=1,3 933 v1 = verl1 * h(i,j) 934 v2 = verl2 * hold(i,j) 935 v3 = verl3 * ( fcell(i,j) - htmp(i,j) ) 936 hnew(i,j) = h(i,j) + ( v1 + v2 + v3 ) * REAL( iforceh(i,j), DP ) 937 ENDDO 938 ENDDO 939 ! 940 END IF 941 942 RETURN 943 END SUBROUTINE cell_verlet 944 945!------------------------------------------------------------------------------! 946 947 subroutine cell_hmove( h, hold, delt, iforceh, fcell ) 948 REAL(DP), intent(out) :: h(3,3) 949 REAL(DP), intent(in) :: hold(3,3), fcell(3,3) 950 REAL(DP), intent(in) :: delt 951 integer, intent(in) :: iforceh(3,3) 952 REAL(DP) :: dt2by2, fac 953 integer :: i, j 954 dt2by2 = 0.5_DP * delt * delt 955 fac = dt2by2 956 do i=1,3 957 do j=1,3 958 h(i,j) = hold(i,j) + fac * iforceh(i,j) * fcell(i,j) 959 end do 960 end do 961 return 962 end subroutine cell_hmove 963 964!------------------------------------------------------------------------------! 965 966 subroutine cell_force( fcell, ainv, stress, omega, press, wmassIN ) 967 USE constants, ONLY : eps8 968 REAL(DP), intent(out) :: fcell(3,3) 969 REAL(DP), intent(in) :: stress(3,3), ainv(3,3) 970 REAL(DP), intent(in) :: omega, press 971 REAL(DP), intent(in), optional :: wmassIN 972 integer :: i, j 973 REAL(DP) :: wmass, fiso 974 IF (.not. present(wmassIN)) THEN 975 wmass = 1.0 976 ELSE 977 wmass = wmassIN 978 END IF 979 do j=1,3 980 do i=1,3 981 fcell(i,j) = ainv(j,1)*stress(i,1) + ainv(j,2)*stress(i,2) + ainv(j,3)*stress(i,3) 982 end do 983 end do 984 do j=1,3 985 do i=1,3 986 fcell(i,j) = fcell(i,j) - ainv(j,i) * press 987 end do 988 end do 989 IF( wmass < eps8 ) & 990 CALL errore( ' movecell ',' cell mass is less than 0 ! ', 1 ) 991 fcell = omega * fcell / wmass 992! added this : 993 IF( isotropic ) THEN 994 ! 995 ! Isotropic force on the cell 996 ! 997 fiso = (fcell(1,1)+fcell(2,2)+fcell(3,3))/3.0_DP 998 do i=1,3 999 fcell(i,i)=fiso 1000 end do 1001 END IF 1002! 1003 return 1004 end subroutine cell_force 1005 1006!------------------------------------------------------------------------------! 1007 1008 subroutine cell_move( hnew, h, hold, delt, iforceh, fcell, frich, tnoseh, vnhh, velh, tsdc ) 1009 REAL(DP), intent(out) :: hnew(3,3) 1010 REAL(DP), intent(in) :: h(3,3), hold(3,3), fcell(3,3) 1011 REAL(DP), intent(in) :: vnhh(3,3), velh(3,3) 1012 integer, intent(in) :: iforceh(3,3) 1013 REAL(DP), intent(in) :: frich, delt 1014 logical, intent(in) :: tnoseh, tsdc 1015 1016 REAL(DP) :: hnos(3,3) 1017 1018 hnew = 0.0 1019 1020 if( tnoseh ) then 1021 hnos = vnhh * velh 1022 else 1023 hnos = 0.0_DP 1024 end if 1025! 1026 IF( tsdc ) THEN 1027 call cell_steepest( hnew, h, delt, iforceh, fcell ) 1028 ELSE 1029 call cell_verlet( hnew, h, hold, delt, iforceh, fcell, frich, tnoseh, hnos ) 1030 END IF 1031 1032 return 1033 end subroutine cell_move 1034 1035!------------------------------------------------------------------------------! 1036 1037 SUBROUTINE cell_gamma( hgamma, ainv, h, velh ) 1038 ! 1039 ! Compute hgamma = g^-1 * dg/dt 1040 ! that enters in the ions equation of motion 1041 ! 1042 IMPLICIT NONE 1043 REAL(DP), INTENT(OUT) :: hgamma(3,3) 1044 REAL(DP), INTENT(IN) :: ainv(3,3), h(3,3), velh(3,3) 1045 REAL(DP) :: gm1(3,3), gdot(3,3) 1046 ! 1047 ! g^-1 inverse of metric tensor = (ht*h)^-1 = ht^-1 * h^-1 1048 ! 1049 gm1 = MATMUL( ainv, TRANSPOSE( ainv ) ) 1050 ! 1051 ! dg/dt = d(ht*h)/dt = dht/dt*h + ht*dh/dt ! derivative of metrix tensor 1052 ! 1053 gdot = MATMUL( TRANSPOSE( velh ), h ) + MATMUL( TRANSPOSE( h ), velh ) 1054 ! 1055 hgamma = MATMUL( gm1, gdot ) 1056 ! 1057 RETURN 1058 END SUBROUTINE cell_gamma 1059 1060!------------------------------------------------------------------------------! 1061 1062 SUBROUTINE cell_update_vel( htp, ht0, htm, delt, velh ) 1063 ! 1064 IMPLICIT NONE 1065 TYPE (boxdimensions) :: htp, ht0, htm 1066 REAL(DP), INTENT(IN) :: delt 1067 REAL(DP), INTENT(OUT) :: velh( 3, 3 ) 1068 1069 velh(:,:) = ( htp%hmat(:,:) - htm%hmat(:,:) ) / ( 2.0d0 * delt ) 1070 htp%gvel = ( htp%g(:,:) - htm%g(:,:) ) / ( 2.0d0 * delt ) 1071 ht0%hvel = velh 1072 1073 RETURN 1074 END SUBROUTINE cell_update_vel 1075 1076 1077!------------------------------------------------------------------------------! 1078 1079 subroutine cell_kinene( ekinh, temphh, velh ) 1080 use constants, only: k_boltzmann_au 1081 implicit none 1082 REAL(DP), intent(out) :: ekinh, temphh(3,3) 1083 REAL(DP), intent(in) :: velh(3,3) 1084 integer :: i,j 1085 ekinh = 0.0_DP 1086 do j=1,3 1087 do i=1,3 1088 ekinh = ekinh + 0.5_DP*wmass*velh(i,j)*velh(i,j) 1089 temphh(i,j) = wmass*velh(i,j)*velh(i,j)/k_boltzmann_au 1090 end do 1091 end do 1092 return 1093 end subroutine cell_kinene 1094 1095!------------------------------------------------------------------------------! 1096 1097 function cell_alat( ) 1098 real(DP) :: cell_alat 1099 if( .NOT. tcell_base_init ) & 1100 call errore( ' cell_alat ', ' alat has not been set ', 1 ) 1101 cell_alat = alat 1102 return 1103 end function cell_alat 1104! 1105!------------------------------------------------------------------------------! 1106 END MODULE cell_base 1107!------------------------------------------------------------------------------! 1108