1! 2! Copyright (C) 2002-2005 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#define __REMOVE_CONSTRAINT_FORCE 9!#define __DEBUG_CONSTRAINTS 10#define __USE_PBC 11! 12!---------------------------------------------------------------------------- 13MODULE constraints_module 14 !---------------------------------------------------------------------------- 15 ! 16 ! ... variables and methods for constraint Molecular Dynamics and 17 ! ... constrained ionic relaxations (the SHAKE algorithm based on 18 ! ... lagrange multipliers) are defined here. 19 ! 20 ! ... most of these variables and methods are also used for meta-dynamics 21 ! ... and free-energy smd : indeed the collective variables are implemented 22 ! ... as constraints. 23 ! 24 ! ... written by Carlo Sbraccia ( 24/02/2004 ) 25 ! 26 ! ... references : 27 ! 28 ! ... 1) M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids, 29 ! ... Clarendon Press - Oxford (1986) 30 ! 31 USE kinds, ONLY : DP 32 USE constants, ONLY : eps32, tpi, fpi 33 USE io_global, ONLY : stdout 34 ! 35 USE basic_algebra_routines 36 ! 37 IMPLICIT NONE 38 ! 39 SAVE 40 ! 41 PRIVATE 42 ! 43 ! ... public methods 44 ! 45 PUBLIC :: init_constraint, & 46 check_constraint, & 47 remove_constr_force, & 48 remove_constr_vec, & 49 deallocate_constraint, & 50 compute_dmax, & 51 pbc, & 52 constraint_grad 53 ! 54 ! 55 ! ... public variables (assigned in the CONSTRAINTS input card) 56 ! 57 PUBLIC :: nconstr, & 58 constr_tol, & 59 constr_type, & 60 constr, & 61 lagrange, & 62 constr_target, & 63 dmax, & 64 gp 65 ! 66 ! ... global variables 67 ! 68 INTEGER :: nconstr=0 69 REAL(DP) :: constr_tol 70 INTEGER, ALLOCATABLE :: constr_type(:) 71 REAL(DP), ALLOCATABLE :: constr(:,:) 72 REAL(DP), ALLOCATABLE :: constr_target(:) 73 REAL(DP), ALLOCATABLE :: lagrange(:) 74 REAL(DP), ALLOCATABLE :: gp(:) 75 REAL(DP) :: dmax 76 ! 77CONTAINS 78 ! 79 ! ... public methods 80 ! 81 !----------------------------------------------------------------------- 82 SUBROUTINE init_constraint( nat, tau, ityp, tau_units ) 83 !----------------------------------------------------------------------- 84 ! 85 ! ... this routine is used to initialize constraints variables and 86 ! ... collective variables (notice that collective variables are 87 ! ... implemented as normal constraints but are read using specific 88 ! ... input variables) 89 ! 90 USE input_parameters, ONLY : nconstr_inp, constr_tol_inp, & 91 constr_type_inp, constr_inp, & 92 constr_target_inp, & 93 constr_target_set, nc_fields 94 ! 95 IMPLICIT NONE 96 ! 97 INTEGER, INTENT(in) :: nat 98 REAL(DP), INTENT(in) :: tau(3,nat) 99 INTEGER, INTENT(in) :: ityp(nat) 100 REAL(DP), INTENT(in) :: tau_units 101 ! 102 INTEGER :: i, j 103 INTEGER :: ia, ia0, ia1, ia2, ia3, n_type_coord1 104 REAL(DP) :: d0(3), d1(3), d2(3) 105 REAL(DP) :: smoothing, r_c 106 INTEGER :: type_coord1, type_coord2 107 REAL(DP) :: dtau(3), norm_dtau 108 REAL(DP) :: k(3), phase, norm_k 109 COMPLEX(DP) :: struc_fac 110 CHARACTER(20),ALLOCATABLE :: tmp_type_inp(:) 111 LOGICAL,ALLOCATABLE :: tmp_target_set(:) 112 REAL(DP),ALLOCATABLE :: tmp_target_inp(:) 113 ! 114 CHARACTER(len=6), EXTERNAL :: int_to_char 115 ! 116 ! 117 nconstr = nconstr_inp 118 constr_tol = constr_tol_inp 119 WRITE(stdout,'(5x,a,i4,a,f12.6)') & 120 'Setting up ',nconstr,' constraints; tolerance:', constr_tol 121 ! 122 ALLOCATE( lagrange( nconstr ) ) 123 ALLOCATE( constr_target( nconstr ) ) 124 ALLOCATE( constr_type( nconstr ) ) 125 ! 126 ALLOCATE( constr( nc_fields, nconstr ) ) 127 ALLOCATE( gp( nconstr ) ) 128 ALLOCATE( tmp_type_inp(nconstr),tmp_target_set(nconstr),tmp_target_inp(nconstr) ) 129 ! 130 ! ... setting constr to 0 to findout which elements have been 131 ! ... set to an atomic index. This is required for CP. 132 ! 133 constr(:,:) = 0.0_DP 134 ! 135 constr(:,1:nconstr) = constr_inp(:,1:nconstr_inp) 136 tmp_type_inp(1:nconstr) = constr_type_inp(1:nconstr_inp) 137 tmp_target_set(1:nconstr) = constr_target_set(1:nconstr_inp) 138 tmp_target_inp(1:nconstr) = constr_target_inp(1:nconstr_inp) 139 ! 140 ! ... set the largest possible distance among two atoms within 141 ! ... the supercell 142 ! 143 IF ( any( tmp_type_inp(:) == 'distance' ) ) CALL compute_dmax() 144 ! 145 ! ... initializations of constr_target values for the constraints : 146 ! 147 DO ia = 1, nconstr 148 ! 149 SELECT CASE ( tmp_type_inp(ia) ) 150 CASE( 'type_coord' ) 151 ! 152 ! ... constraint on global coordination-number, i.e. the average 153 ! ... number of atoms of type B surrounding the atoms of type A 154 ! 155 constr_type(ia) = 1 156 IF ( tmp_target_set(ia) ) THEN 157 constr_target(ia) = tmp_target_inp(ia) 158 ELSE 159 CALL set_type_coord( ia ) 160 ENDIF 161 ! 162 WRITE(stdout,'(7x,i3,a,i3,a,i2,a,2f12.6,a,f12.6)') & 163 ia,') type #',int(constr_inp(1,ia)) ,' coordination wrt type:', int(constr(2,ia)), & 164 ' cutoff distance and smoothing:', constr(3:4,ia), & 165 '; target:', constr_target(ia) 166 ! 167 CASE( 'atom_coord' ) 168 ! 169 ! ... constraint on local coordination-number, i.e. the average 170 ! ... number of atoms of type A surrounding a specific atom 171 ! 172 constr_type(ia) = 2 173 IF ( tmp_target_set(ia) ) THEN 174 constr_target(ia) = tmp_target_inp(ia) 175 ELSE 176 CALL set_atom_coord( ia ) 177 ENDIF 178 ! 179 WRITE(stdout,'(7x,i3,a,i3,a,i2,a,2f12.6,a,f12.6)') & 180 ia,') atom #',int(constr_inp(1,ia)) ,' coordination wrt type:', int(constr(2,ia)), & 181 ' cutoff distance and smoothing:', constr(3:4,ia), & 182 '; target:', constr_target(ia) 183 ! 184 CASE( 'distance' ) 185 ! 186 constr_type(ia) = 3 187 IF ( tmp_target_set(ia) ) THEN 188 constr_target(ia) = tmp_target_inp(ia) 189 ELSE 190 CALL set_distance( ia ) 191 ENDIF 192 ! 193 IF ( constr_target(ia) > dmax ) THEN 194 ! 195 WRITE( stdout, '(/,5X,"target = ",F12.8,/, & 196 & 5X,"dmax = ",F12.8)' ) & 197 constr_target(ia), dmax 198 CALL errore( 'init_constraint', 'the target for constraint ' //& 199 & trim( int_to_char( ia ) ) // ' is larger than ' //& 200 & 'the largest possible value', 1 ) 201 ! 202 ENDIF 203 ! 204 WRITE(stdout,'(7x,i3,a,2i3,a,f12.6)') & 205 ia,') distance between atoms: ', int(constr(1:2,ia)), '; target:', constr_target(ia) 206 ! 207 CASE( 'planar_angle' ) 208 ! 209 ! ... constraint on planar angle (for the notation used here see 210 ! ... Appendix C of the Allen-Tildesley book) 211 ! 212 constr_type(ia) = 4 213 IF ( tmp_target_set(ia) ) THEN 214 ! 215 ! ... the input value of target for the torsional angle (given 216 ! ... in degrees) is converted to the cosine of the angle 217 ! 218 constr_target(ia) = tmp_target_inp(ia) 219 ELSE 220 CALL set_planar_angle( ia ) 221 ENDIF 222 ! 223 WRITE(stdout, '(7x,i3,a,3i3,a,f12.6)') & 224 ia,') planar angle between atoms: ', int(constr(1:3,ia)), '; target:', constr_target(ia) 225 ! 226 CASE( 'torsional_angle' ) 227 ! 228 ! ... constraint on torsional angle (for the notation used here 229 ! ... see Appendix C of the Allen-Tildesley book) 230 ! 231 constr_type(ia) = 5 232 IF ( tmp_target_set(ia) ) THEN 233 ! 234 ! ... the input value of target for the torsional angle (given 235 ! ... in degrees) is converted to the cosine of the angle 236 ! 237 constr_target(ia) = tmp_target_inp(ia) 238 ELSE 239 CALL set_torsional_angle( ia ) 240 ENDIF 241 ! 242 WRITE(stdout, '(7x,i3,a,4i3,a,f12.6)') & 243 ia,') torsional angle between atoms: ', int(constr(1:4,ia)), '; target:', constr_target(ia) 244 ! 245 CASE( 'struct_fac' ) 246 ! 247 ! ... constraint on structure factor at a given k-vector 248 ! 249 constr_type(ia) = 6 250 IF ( tmp_target_set(ia) ) THEN 251 constr_target(ia) = tmp_target_inp(ia) 252 ELSE 253 CALL set_structure_factor( ia ) 254 ENDIF 255 ! 256 CASE( 'sph_struct_fac' ) 257 ! 258 ! ... constraint on spherical average of the structure factor for 259 ! ... a given k-vector of norm k 260 ! 261 constr_type(ia) = 7 262 IF ( tmp_target_set(ia) ) THEN 263 constr_target(ia) = tmp_target_inp(ia) 264 ELSE 265 CALL set_sph_structure_factor( ia ) 266 ENDIF 267 ! 268 CASE( 'bennett_proj' ) 269 ! 270 ! ... constraint on the projection onto a given direction of the 271 ! ... vector defined by the position of one atom minus the center 272 ! ... of mass of the others 273 ! ... ( Ch.H. Bennett in Diffusion in Solids, Recent Developments, 274 ! ... Ed. by A.S. Nowick and J.J. Burton, New York 1975 ) 275 ! 276 constr_type(ia) = 8 277 IF ( tmp_target_set(ia) ) THEN 278 constr_target(ia) = tmp_target_inp(ia) 279 ELSE 280 CALL set_bennett_proj( ia ) 281 ENDIF 282 ! 283 CASE DEFAULT 284 ! 285 CALL errore( 'init_constraint', & 286 'collective-variable or constrait type not implemented', 1 ) 287 ! 288 END SELECT 289 ! 290 ENDDO 291 ! 292 DEALLOCATE( tmp_type_inp,tmp_target_set,tmp_target_inp ) 293 ! 294 RETURN 295 ! 296 CONTAINS 297 ! 298 !------------------------------------------------------------------- 299 SUBROUTINE set_type_coord( ia ) 300 !------------------------------------------------------------------- 301 ! 302 INTEGER, INTENT(in) :: ia 303 ! 304 type_coord1 = anint( constr(1,ia) ) 305 type_coord2 = anint( constr(2,ia) ) 306 ! 307 r_c = constr(3,ia) 308 ! 309 smoothing = 1.0_DP / constr(4,ia) 310 ! 311 constr_target(ia) = 0.0_DP 312 ! 313 n_type_coord1 = 0 314 ! 315 DO ia1 = 1, nat 316 ! 317 IF ( ityp(ia1) /= type_coord1 ) CYCLE 318 ! 319 DO ia2 = 1, nat 320 ! 321 IF ( ia2 == ia1 ) CYCLE 322 ! 323 IF ( ityp(ia2) /= type_coord2 ) CYCLE 324 ! 325 dtau(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 326 ! 327 norm_dtau = norm( dtau(:) ) 328 ! 329 constr_target(ia) = constr_target(ia) + 1.0_DP / & 330 ( exp( smoothing*( norm_dtau - r_c ) ) + 1.0_DP ) 331 ! 332 ENDDO 333 ! 334 n_type_coord1 = n_type_coord1 + 1 335 ! 336 ENDDO 337 ! 338 constr_target(ia) = constr_target(ia) / dble( n_type_coord1 ) 339 ! 340 END SUBROUTINE set_type_coord 341 ! 342 !------------------------------------------------------------------- 343 SUBROUTINE set_atom_coord( ia ) 344 !------------------------------------------------------------------- 345 ! 346 INTEGER, INTENT(in) :: ia 347 ! 348 ia1 = anint( constr(1,ia) ) 349 type_coord1 = anint( constr(2,ia) ) 350 ! 351 r_c = constr(3,ia) 352 ! 353 smoothing = 1.0_DP / constr(4,ia) 354 ! 355 constr_target(ia) = 0.0_DP 356 ! 357 DO ia2 = 1, nat 358 ! 359 IF ( ia2 == ia1 ) CYCLE 360 ! 361 IF ( ityp(ia2) /= type_coord1 ) CYCLE 362 ! 363 dtau(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 364 ! 365 norm_dtau = norm( dtau(:) ) 366 ! 367 constr_target(ia) = constr_target(ia) + 1.0_DP / & 368 ( exp( smoothing*( norm_dtau - r_c ) ) + 1.0_DP ) 369 ! 370 ENDDO 371 ! 372 END SUBROUTINE set_atom_coord 373 ! 374 !------------------------------------------------------------------- 375 SUBROUTINE set_distance( ia ) 376 !------------------------------------------------------------------- 377 ! 378 INTEGER, INTENT(in) :: ia 379 ! 380 ia1 = anint( constr(1,ia) ) 381 ia2 = anint( constr(2,ia) ) 382 ! 383 dtau(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 384 ! 385 constr_target(ia) = norm( dtau(:) ) 386 ! 387 END SUBROUTINE set_distance 388 ! 389 !------------------------------------------------------------------- 390 SUBROUTINE set_planar_angle( ia ) 391 !------------------------------------------------------------------- 392 ! 393 INTEGER, INTENT(in) :: ia 394 ! 395 ia0 = anint( constr(1,ia) ) 396 ia1 = anint( constr(2,ia) ) 397 ia2 = anint( constr(3,ia) ) 398 ! 399 d0(:) = pbc( ( tau(:,ia0) - tau(:,ia1) )*tau_units ) 400 d1(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 401 ! 402 d0(:) = d0(:) / norm( d0(:) ) 403 d1(:) = d1(:) / norm( d1(:) ) 404 ! 405 constr_target(ia) = acos(- d0(:) .dot. d1(:))*360.0_DP/tpi 406 ! 407 END SUBROUTINE set_planar_angle 408 ! 409 !------------------------------------------------------------------- 410 SUBROUTINE set_torsional_angle( ia ) 411 !------------------------------------------------------------------- 412 ! 413 INTEGER, INTENT(in) :: ia 414 REAL(DP) :: x01(3),x12(3),phi 415 ! 416 ia0 = anint( constr(1,ia) ) 417 ia1 = anint( constr(2,ia) ) 418 ia2 = anint( constr(3,ia) ) 419 ia3 = anint( constr(4,ia) ) 420 ! 421 d0(:) = pbc( ( tau(:,ia0) - tau(:,ia1) )*tau_units ) 422 d1(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 423 d2(:) = pbc( ( tau(:,ia2) - tau(:,ia3) )*tau_units ) 424 ! 425 x01(:) = cross(d0,d1) 426 x12(:) = cross(d1,d2) 427 ! 428 IF((x01.dot.x01)<eps32 .or. (x12.dot.x12)<eps32)THEN 429 write(stdout,*)'torsional angle constraint #',ia,' contains collinear atoms' 430 CALL errore('set_torsional_angle','collinear atoms in torsional angle constraint', 1) 431 ENDIF 432 ! 433 phi = atan2(sqrt(d1.dot.d1)*d0.dot.x12 , x01.dot.x12) 434 ! 435 constr_target(ia) = phi*360.0_DP/tpi 436 ! 437 END SUBROUTINE set_torsional_angle 438 ! 439 !------------------------------------------------------------------- 440 SUBROUTINE set_structure_factor( ia ) 441 !------------------------------------------------------------------- 442 ! 443 INTEGER, INTENT(in) :: ia 444 ! 445 k(1) = constr(1,ia) * tpi / tau_units 446 k(2) = constr(2,ia) * tpi / tau_units 447 k(3) = constr(3,ia) * tpi / tau_units 448 ! 449 struc_fac = ( 0.0_DP, 0.0_DP ) 450 ! 451 DO i = 1, nat 452 ! 453 dtau(:) = pbc( ( tau(:,i) - tau(:,1) )*tau_units ) 454 ! 455 phase = k(:) .dot. dtau(:) 456 ! 457 struc_fac = struc_fac + cmplx( cos(phase), sin(phase), kind=DP ) 458 ! 459 ENDDO 460 ! 461 constr_target(ia) = conjg( struc_fac )*struc_fac / dble( nat*nat ) 462 ! 463 END SUBROUTINE set_structure_factor 464 ! 465 !------------------------------------------------------------------- 466 SUBROUTINE set_sph_structure_factor( ia ) 467 !------------------------------------------------------------------- 468 ! 469 INTEGER, INTENT(in) :: ia 470 ! 471 norm_k = constr(1,ia)*tpi/tau_units 472 ! 473 constr_target(ia) = 0.0_DP 474 ! 475 DO i = 1, nat - 1 476 ! 477 DO j = i + 1, nat 478 ! 479 dtau(:) = pbc( ( tau(:,i) - tau(:,j) )*tau_units ) 480 ! 481 norm_dtau = norm( dtau(:) ) 482 ! 483 phase = norm_k*norm_dtau 484 ! 485 IF ( phase < eps32 ) THEN 486 ! 487 constr_target(ia) = constr_target(ia) + 1.0_DP 488 ! 489 ELSE 490 ! 491 constr_target(ia) = constr_target(ia) + sin( phase ) / phase 492 ! 493 ENDIF 494 ! 495 ENDDO 496 ! 497 ENDDO 498 ! 499 constr_target(ia) = 2.0_DP * fpi * constr_target(ia) / dble( nat ) 500 ! 501 END SUBROUTINE set_sph_structure_factor 502 ! 503 !------------------------------------------------------------------- 504 SUBROUTINE set_bennett_proj( ia ) 505 !------------------------------------------------------------------- 506 ! 507 INTEGER, INTENT(in) :: ia 508 ! 509 ia0 = anint( constr(1,ia) ) 510 ! 511 d0(:) = tau(:,ia0) 512 d1(:) = sum( tau(:,:), dim = 2 ) 513 ! 514 d1(:) = pbc( ( d1(:) - d0(:) )*tau_units ) / dble( nat - 1 ) - & 515 pbc( d0(:)*tau_units ) 516 ! 517 d2(:) = constr(2:4,ia) 518 ! 519 constr_target(ia) = ( d1(:) .dot. d2(:) ) / tau_units 520 ! 521 END SUBROUTINE set_bennett_proj 522 ! 523 END SUBROUTINE init_constraint 524 ! 525 !----------------------------------------------------------------------- 526 SUBROUTINE constraint_grad( idx, nat, tau, & 527 if_pos, ityp, tau_units, g, dg ) 528 !----------------------------------------------------------------------- 529 ! 530 ! ... this routine computes the value of the constraint equation and 531 ! ... the corresponding constraint gradient 532 ! 533 IMPLICIT NONE 534 ! 535 INTEGER, INTENT(in) :: idx 536 INTEGER, INTENT(in) :: nat 537 REAL(DP), INTENT(in) :: tau(:,:) 538 INTEGER, INTENT(in) :: if_pos(:,:) 539 INTEGER, INTENT(in) :: ityp(:) 540 REAL(DP), INTENT(in) :: tau_units 541 REAL(DP), INTENT(out) :: dg(:,:) 542 REAL(DP), INTENT(out) :: g 543 ! 544 INTEGER :: i, j 545 INTEGER :: ia, ia0, ia1, ia2, ia3, n_type_coord1 546 REAL(DP) :: d0(3), d1(3), d2(3), n1, x01(3), x12(3), x20(3), phi, X 547 REAL(DP) :: s012,x01x12,d0phi(3),d1phi(3),d2phi(3) 548 REAL(DP) :: inv_den 549 REAL(DP) :: C00, C01, C11 550 REAL(DP) :: smoothing, r_c 551 INTEGER :: type_coord1, type_coord2 552 REAL(DP) :: dtau(3), norm_dtau, norm_dtau_sq, expo 553 REAL(DP) :: r0(3), r1(3), r2(3), ri(3), k(3), phase, ksin(3), norm_k, sinxx 554 COMPLEX(DP) :: struc_fac 555 ! 556 dg(:,:) = 0.0_DP 557 ! 558 SELECT CASE ( constr_type(idx) ) 559 CASE( 1 ) 560 ! 561 ! ... constraint on global coordination 562 ! 563 type_coord1 = anint( constr(1,idx) ) 564 type_coord2 = anint( constr(2,idx) ) 565 ! 566 r_c = constr(3,idx) 567 ! 568 smoothing = 1.0_DP / constr(4,idx) 569 ! 570 g = 0.0_DP 571 ! 572 n_type_coord1 = 0 573 ! 574 DO ia1 = 1, nat 575 ! 576 IF ( ityp(ia1) /= type_coord1 ) CYCLE 577 ! 578 DO ia2 = 1, nat 579 ! 580 IF ( ia2 == ia1 ) CYCLE 581 ! 582 IF ( ityp(ia2) /= type_coord2 ) CYCLE 583 ! 584 dtau(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 585 ! 586 norm_dtau = norm( dtau(:) ) 587 ! 588 dtau(:) = dtau(:) / norm_dtau 589 ! 590 expo = exp( smoothing*( norm_dtau - r_c ) ) 591 ! 592 g = g + 1.0_DP / ( expo + 1.0_DP ) 593 ! 594 dtau(:) = dtau(:) * smoothing*expo / ( expo + 1.0_DP )**2 595 ! 596 dg(:,ia2) = dg(:,ia2) + dtau(:) 597 dg(:,ia1) = dg(:,ia1) - dtau(:) 598 ! 599 ENDDO 600 ! 601 n_type_coord1 = n_type_coord1 + 1 602 ! 603 ENDDO 604 ! 605 g = g / dble( n_type_coord1 ) 606 dg = dg / dble( n_type_coord1 ) 607 ! 608 g = ( g - constr_target(idx) ) 609 ! 610 CASE( 2 ) 611 ! 612 ! ... constraint on local coordination 613 ! 614 ia = anint( constr(1,idx) ) 615 type_coord1 = anint( constr(2,idx) ) 616 ! 617 r_c = constr(3,idx) 618 ! 619 smoothing = 1.0_DP / constr(4,idx) 620 ! 621 g = 0.0_DP 622 ! 623 DO ia1 = 1, nat 624 ! 625 IF ( ia1 == ia ) CYCLE 626 ! 627 IF ( ityp(ia1) /= type_coord1 ) CYCLE 628 ! 629 dtau(:) = pbc( ( tau(:,ia) - tau(:,ia1) )*tau_units ) 630 ! 631 norm_dtau = norm( dtau(:) ) 632 ! 633 dtau(:) = dtau(:) / norm_dtau 634 ! 635 expo = exp( smoothing*( norm_dtau - r_c ) ) 636 ! 637 g = g + 1.0_DP / ( expo + 1.0_DP ) 638 ! 639 dtau(:) = dtau(:) * smoothing * expo / ( expo + 1.0_DP )**2 640 ! 641 dg(:,ia1) = dg(:,ia1) + dtau(:) 642 dg(:,ia) = dg(:,ia) - dtau(:) 643 ! 644 ENDDO 645 ! 646 g = ( g - constr_target(idx) ) 647 ! 648 CASE( 3 ) 649 ! 650 ! ... constraint on distances 651 ! 652 ia1 = anint( constr(1,idx) ) 653 ia2 = anint( constr(2,idx) ) 654 ! 655 dtau(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 656 ! 657 norm_dtau = norm( dtau(:) ) 658 ! 659 g = ( norm_dtau - constr_target(idx) ) 660 ! 661 dg(:,ia1) = dtau(:) / norm_dtau 662 ! 663 dg(:,ia2) = - dg(:,ia1) 664 ! 665 CASE( 4 ) 666 ! 667 ! ... constraint on planar angles (for the notation used here see 668 ! ... Appendix C of the Allen-Tildesley book) 669 ! 670 ia0 = anint( constr(1,idx) ) 671 ia1 = anint( constr(2,idx) ) 672 ia2 = anint( constr(3,idx) ) 673 ! 674 d0(:) = pbc( ( tau(:,ia0) - tau(:,ia1) )*tau_units ) 675 d1(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 676 ! 677 C00 = d0(:) .dot. d0(:) 678 C01 = d0(:) .dot. d1(:) 679 C11 = d1(:) .dot. d1(:) 680 ! 681 inv_den = 1.0_DP / sqrt( C00*C11 ) 682 ! 683 g = ( acos(- C01 * inv_den)*360.d0/tpi - constr_target(idx) ) 684 ! 685 ! d/dx acos(x) = -1/sqrt(1-x**2) 686 ! d/dx acos(-x)*360/tpi = (360/tpi)/sqrt(1-x**2)) 687 ! 688 X = (360.d0/tpi)/sqrt(1-(C01 * inv_den)**2) 689 dg(:,ia0) = X * ( d1(:) - C01/C00*d0(:) ) * inv_den 690 dg(:,ia2) = X * ( C01/C11*d1(:) - d0(:) ) * inv_den 691 dg(:,ia1) = - dg(:,ia0) - dg(:,ia2) 692 ! 693 CASE( 5 ) 694 ! 695 ! ... constraint on torsional angle (for the notation used here 696 ! ... see Appendix C of the Allen-Tildesley book) 697 ! 698 ia0 = anint( constr(1,idx) ) 699 ia1 = anint( constr(2,idx) ) 700 ia2 = anint( constr(3,idx) ) 701 ia3 = anint( constr(4,idx) ) 702 ! 703 r0(:) = pbc( ( tau(:,ia0) - tau(:,ia1) )*tau_units ) 704 r1(:) = pbc( ( tau(:,ia1) - tau(:,ia2) )*tau_units ) 705 r2(:) = pbc( ( tau(:,ia2) - tau(:,ia3) )*tau_units ) 706 n1 = sqrt(r1.dot.r1) 707 ! 708 x01(:) = cross(r0,r1) 709 x12(:) = cross(r1,r2) 710 x20(:) = cross(r2,r0) 711 ! 712 s012 = r0.dot.x12 713 x01x12 = x01.dot.x12 714 ! 715 phi = atan2(n1*s012 , x01x12) 716 ! 717 g = phi*360.0_DP/tpi - constr_target(idx) 718 g = modulo(g+180.0_DP,360.0_DP)-180.0_DP 719 ! 720 ! d/dx atan(x) = 1/1+x**2 721 ! d/dy atan2(y,x) = x/(x**2+y**2) 722 ! d/dx atan2(y,x) = -y/(x**2+y**2) 723 ! d(atan2(A,B)) = (BdA-AdB)/(A**2+B**2) 724 ! 725 ! dd0(:,ia0) = 1 ; dd1(:,ia0) = 0 ; dd2(:,ia0) = 0 726 ! dd0(:,ia1) = -1 ; dd1(:,ia1) = 1 ; dd2(:,ia1) = 0 727 ! dd0(:,ia2) = 0 ; dd1(:,ia2) = -1 ; dd2(:,ia2) = 1 728 ! dd0(:,ia3) = 0 ; dd1(:,ia3) = 0 ; dd2(:,ia3) = -1 729 ! 730 ! d(s012) / d(r0) = x12 731 ! d(s012) / d(r1) = x20 732 ! d(s012) / d(r2) = x01 733 ! 734 ! d(x01x12) / d(r0) = r1 x x12 735 ! d(x01x12) / d(r1) = r2 x x01 + x12 x r0 736 ! d(x01x12) / d(r2) = x01 x r1 737 ! 738 ! d(n1) / d(r1) = r1 / n1 739 ! 740 ! d(phi) = (x01x12 * d(n1 * s012) - n1*s012 * d(x01x12) 741 ! / 742 ! (x01x12 ** 2 + n1*s012 ** 2) 743 ! 744 ! d(phi)/d(d0) = (x01x12 * n1*x12 - n1*s012 * r1 x x12) 745 ! / DENOM 746 ! d(phi)/d(d1) = (x01x12 * (n1*x20 + d1/n1 * s012) - n1*s012 * (r2 x x01 + x12 x r0)) 747 ! / DENOM 748 ! d(phi)/d(d2) = (x01x12 * n1*x01 - n1*s012 * x01 x r1) 749 ! / DENOM 750 ! 751 inv_den = 1.0_DP / (x01x12 ** 2 + n1*s012 ** 2) 752 d0phi(:) = (x01x12 * n1*x12 - n1*s012 * cross(r1,x12)) * inv_den 753 d1phi(:) = (x01x12 * (n1*x20 + d1/n1 * s012) - n1*s012 * (cross(r2,x01) + cross(x12,r0))) * inv_den 754 d2phi(:) = (x01x12 * n1*x01 - n1*s012 * cross(x01,r1)) * inv_den 755 ! 756 dg(:,ia0) = d0phi*360.0_DP/tpi 757 dg(:,ia1) = (d1phi-d0phi)*360.0_DP/tpi 758 dg(:,ia2) = (d2phi-d1phi)*360.0_DP/tpi 759 dg(:,ia3) = (-d2phi)*360.0_DP/tpi 760 ! 761 CASE( 6 ) 762 ! 763 ! ... constraint on structure factor at a given k vector 764 ! 765 k(1) = constr(1,idx)*tpi/tau_units 766 k(2) = constr(2,idx)*tpi/tau_units 767 k(3) = constr(3,idx)*tpi/tau_units 768 ! 769 struc_fac = ( 1.0_DP, 0.0_DP ) 770 ! 771 r0(:) = tau(:,1) 772 ! 773 DO i = 1, nat - 1 774 ! 775 dtau(:) = pbc( ( tau(:,i+1) - r0(:) )*tau_units ) 776 ! 777 phase = k(1)*dtau(1) + k(2)*dtau(2) + k(3)*dtau(3) 778 ! 779 struc_fac = struc_fac + cmplx( cos(phase), sin(phase), kind=DP ) 780 ! 781 ri(:) = tau(:,i) 782 ! 783 DO j = i + 1, nat 784 ! 785 dtau(:) = pbc( ( tau(:,j) - ri(:) )*tau_units ) 786 ! 787 phase = k(1)*dtau(1) + k(2)*dtau(2) + k(3)*dtau(3) 788 ! 789 ksin(:) = k(:)*sin( phase ) 790 ! 791 dg(:,i) = dg(:,i) + ksin(:) 792 dg(:,j) = dg(:,j) - ksin(:) 793 ! 794 ENDDO 795 ! 796 ENDDO 797 ! 798 g = ( conjg( struc_fac )*struc_fac ) / dble( nat*nat ) 799 ! 800 g = ( g - constr_target(idx) ) 801 ! 802 dg(:,:) = dg(:,:)*2.0_DP/dble( nat*nat ) 803 ! 804 CASE( 7 ) 805 ! 806 ! ... constraint on spherical average of the structure factor for 807 ! ... a given k-vector of norm k 808 ! 809 norm_k = constr(1,idx)*tpi/tau_units 810 ! 811 g = 0.0_DP 812 ! 813 DO i = 1, nat - 1 814 ! 815 ri(:) = tau(:,i) 816 ! 817 DO j = i + 1, nat 818 ! 819 dtau(:) = pbc( ( ri(:) - tau(:,j) )*tau_units ) 820 ! 821 norm_dtau_sq = dtau(1)**2 + dtau(2)**2 + dtau(3)**2 822 ! 823 norm_dtau = sqrt( norm_dtau_sq ) 824 ! 825 phase = norm_k * norm_dtau 826 ! 827 IF ( phase < eps32 ) THEN 828 ! 829 g = g + 1.0_DP 830 ! 831 ELSE 832 ! 833 sinxx = sin( phase ) / phase 834 ! 835 g = g + sinxx 836 ! 837 dtau(:) = dtau(:) / norm_dtau_sq*( cos( phase ) - sinxx ) 838 ! 839 dg(:,i) = dg(:,i) + dtau(:) 840 dg(:,j) = dg(:,j) - dtau(:) 841 ! 842 ENDIF 843 ! 844 ENDDO 845 ! 846 ENDDO 847 ! 848 g = ( 2.0_DP*fpi*g / dble( nat ) - constr_target(idx) ) 849 ! 850 dg(:,:) = 4.0_DP*fpi*dg(:,:) / dble( nat ) 851 ! 852 CASE( 8 ) 853 ! 854 ! ... constraint on Bennett projection 855 ! 856 ia0 = anint( constr(1,idx) ) 857 ! 858 d0(:) = tau(:,ia0) 859 d1(:) = sum( tau(:,:), dim = 2 ) 860 ! 861 d1(:) = pbc( ( d1(:) - d0(:) )*tau_units ) / dble( nat - 1 ) - & 862 pbc( d0(:)*tau_units ) 863 ! 864 d2(:) = constr(2:4,idx) 865 ! 866 g = ( d1(:) .dot. d2(:) ) / tau_units - constr_target( idx ) 867 ! 868 dg = 0.0_DP 869 ! 870 C00 = ( 1.0_DP / dble( nat - 1 ) ) / tau_units 871 C01 = -1.0_DP / tau_units 872 ! 873 DO i = 1, nat 874 ! 875 dg(:,i) = d2(:)*C00 876 ! 877 ENDDO 878 ! 879 dg(:,ia0) = d2(:)*C01 880 ! 881 END SELECT 882 ! 883 dg(:,:) = dg(:,:)*dble( if_pos(:,:) ) 884 ! 885 RETURN 886 ! 887 END SUBROUTINE constraint_grad 888 ! 889 !----------------------------------------------------------------------- 890 SUBROUTINE check_constraint( nat, taup, tau0, & 891 force, if_pos, ityp, tau_units, dt, massconv ) 892 !----------------------------------------------------------------------- 893 ! 894 ! ... update taup (predicted positions) so that the constraint equation 895 ! ... g=0 is satisfied, using the recursion formula: 896 ! 897 ! ... g(taup) 898 ! ... taup = taup - ----------------------- * dg(tau0) 899 ! ... M^-1<dg(taup)|dg(tau0)> 900 ! 901 ! ... in normal cases the constraint equation should be satisfied at 902 ! ... the very first iteration. 903 ! 904 USE ions_base, ONLY : amass 905 ! 906 IMPLICIT NONE 907 ! 908 INTEGER, INTENT(in) :: nat 909 REAL(DP), INTENT(inout) :: taup(3,nat) 910 REAL(DP), INTENT(in) :: tau0(3,nat) 911 INTEGER, INTENT(in) :: if_pos(3,nat) 912 REAL(DP), INTENT(inout) :: force(3,nat) 913 INTEGER, INTENT(in) :: ityp(nat) 914 REAL(DP), INTENT(in) :: tau_units 915 REAL(DP), INTENT(in) :: dt 916 REAL(DP), INTENT(in) :: massconv 917 ! 918 INTEGER :: na, i, idx, dim 919 REAL(DP), ALLOCATABLE :: dgp(:,:), dg0(:,:,:) 920 REAL(DP) :: g0 921 REAL(DP) :: lambda, fac, invdtsq 922 LOGICAL, ALLOCATABLE :: ltest(:) 923 LOGICAL :: global_test 924 INTEGER, PARAMETER :: maxiter = 100 925 ! 926 REAL(DP), EXTERNAL :: ddot 927 ! 928 ! 929 ALLOCATE( dgp( 3, nat ) ) 930 ALLOCATE( dg0( 3, nat, nconstr ) ) 931 ! 932 ALLOCATE( ltest( nconstr ) ) 933 ! 934 invdtsq = 1.0_DP / dt**2 935 ! 936 dim = 3*nat 937 ! 938 DO idx = 1, nconstr 939 ! 940 CALL constraint_grad( idx, nat, tau0, & 941 if_pos, ityp, tau_units, g0, dg0(:,:,idx) ) 942 ! 943 ENDDO 944 ! 945 outer_loop: DO i = 1, maxiter 946 ! 947 inner_loop: DO idx = 1, nconstr 948 ! 949 ltest(idx) = .false. 950 ! 951 CALL constraint_grad( idx, nat, taup, & 952 if_pos, ityp, tau_units, gp(idx), dgp ) 953 ! 954 ! ... check if gp = 0 955 ! 956#if defined (__DEBUG_CONSTRAINTS) 957 WRITE( stdout, '(2(2X,I3),F12.8)' ) i, idx, abs( gp(idx) ) 958#endif 959 ! 960 IF ( abs( gp(idx) ) < constr_tol ) THEN 961 ! 962 ltest(idx) = .true. 963 ! 964 CYCLE inner_loop 965 ! 966 ENDIF 967 ! 968 ! ... if gp <> 0 find new taup and check again 969 ! ... ( gp is in bohr and taup in tau_units ) 970 ! 971 DO na = 1, nat 972 ! 973 dgp(:,na) = dgp(:,na) / ( amass(ityp(na))*massconv ) 974 ! 975 ENDDO 976 ! 977 lambda = gp(idx) / ddot( dim, dgp, 1, dg0(:,:,idx), 1 ) 978 ! 979 DO na = 1, nat 980 ! 981 fac = amass(ityp(na))*massconv*tau_units 982 ! 983 taup(:,na) = taup(:,na) - lambda*dg0(:,na,idx)/fac 984 ! 985 ENDDO 986 ! 987 lagrange(idx) = lagrange(idx) + lambda*invdtsq 988 ! 989 force(:,:) = force(:,:) - lambda*dg0(:,:,idx)*invdtsq 990 ! 991 ENDDO inner_loop 992 ! 993 global_test = all( ltest(:) ) 994 ! 995 ! ... all constraints are satisfied 996 ! 997 IF ( global_test ) exit outer_loop 998 ! 999 ENDDO outer_loop 1000 ! 1001 IF ( .not. global_test ) THEN 1002 ! 1003 ! ... error messages 1004 ! 1005 WRITE( stdout, '(/,5X,"Number of step(s): ",I3)') min( i, maxiter ) 1006 WRITE( stdout, '(/,5X,"constr_target convergence: ")' ) 1007 ! 1008 DO i = 1, nconstr 1009 ! 1010 WRITE( stdout, '(5X,"constr # ",I3,2X,L1,3(2X,F16.10))' ) & 1011 i, ltest(i), abs( gp(i) ), constr_tol, constr_target(i) 1012 ! 1013 ENDDO 1014 ! 1015 CALL errore( 'check_constraint', & 1016 'on some constraint g = 0 is not satisfied', 1 ) 1017 ! 1018 ENDIF 1019 ! 1020 DEALLOCATE( dgp ) 1021 DEALLOCATE( dg0 ) 1022 DEALLOCATE( ltest ) 1023 ! 1024 RETURN 1025 ! 1026 END SUBROUTINE check_constraint 1027 ! 1028 !----------------------------------------------------------------------- 1029 SUBROUTINE remove_constr_force( nat, tau, & 1030 if_pos, ityp, tau_units, force ) 1031 !----------------------------------------------------------------------- 1032 ! 1033 ! ... the component of the force that is orthogonal to the 1034 ! ... ipersurface defined by the constraint equations is removed 1035 ! ... and the corresponding value of the lagrange multiplier computed 1036 ! 1037 IMPLICIT NONE 1038 ! 1039 INTEGER, INTENT(in) :: nat 1040 REAL(DP), INTENT(in) :: tau(:,:) 1041 INTEGER, INTENT(in) :: if_pos(:,:) 1042 INTEGER, INTENT(in) :: ityp(:) 1043 REAL(DP), INTENT(in) :: tau_units 1044 REAL(DP), INTENT(inout) :: force(:,:) 1045 ! 1046 INTEGER :: i, j, dim 1047 REAL(DP) :: g, ndg, dgidgj 1048 REAL(DP) :: norm_before, norm_after 1049 REAL(DP), ALLOCATABLE :: dg(:,:,:) 1050 REAL(DP), ALLOCATABLE :: dg_matrix(:,:) 1051 INTEGER, ALLOCATABLE :: iwork(:) 1052 ! 1053 REAL(DP), EXTERNAL :: ddot, dnrm2 1054 ! 1055 ! 1056 dim = 3*nat 1057 ! 1058 lagrange(:) = 0.0_DP 1059 ! 1060#if defined (__REMOVE_CONSTRAINT_FORCE) 1061 ! 1062 norm_before = dnrm2( 3*nat, force, 1 ) 1063 ! 1064 ALLOCATE( dg( 3, nat, nconstr ) ) 1065 ! 1066 IF ( nconstr == 1 ) THEN 1067 ! 1068 CALL constraint_grad( 1, nat, tau, & 1069 if_pos, ityp, tau_units, g, dg(:,:,1) ) 1070 ! 1071 lagrange(1) = ddot( dim, force, 1, dg(:,:,1), 1 ) 1072 ! 1073 ndg = ddot( dim, dg(:,:,1), 1, dg(:,:,1), 1 ) 1074 ! 1075 force(:,:) = force(:,:) - lagrange(1)*dg(:,:,1)/ndg 1076 ! 1077 ELSE 1078 ! 1079 ALLOCATE( dg_matrix( nconstr, nconstr ) ) 1080 ALLOCATE( iwork( nconstr ) ) 1081 ! 1082 DO i = 1, nconstr 1083 ! 1084 CALL constraint_grad( i, nat, tau, & 1085 if_pos, ityp, tau_units, g, dg(:,:,i) ) 1086 ! 1087 ENDDO 1088 ! 1089 DO i = 1, nconstr 1090 ! 1091 dg_matrix(i,i) = ddot( dim, dg(:,:,i), 1, dg(:,:,i), 1 ) 1092 ! 1093 lagrange(i) = ddot( dim, force, 1, dg(:,:,i), 1 ) 1094 ! 1095 DO j = i + 1, nconstr 1096 ! 1097 dgidgj = ddot( dim, dg(:,:,i), 1, dg(:,:,j), 1 ) 1098 ! 1099 dg_matrix(i,j) = dgidgj 1100 dg_matrix(j,i) = dgidgj 1101 ! 1102 ENDDO 1103 ! 1104 ENDDO 1105 ! 1106 CALL DGESV( nconstr, 1, dg_matrix, & 1107 nconstr, iwork, lagrange, nconstr, i ) 1108 ! 1109 IF ( i /= 0 ) & 1110 CALL errore( 'remove_constr_force', & 1111 'error in the solution of the linear system', i ) 1112 ! 1113 DO i = 1, nconstr 1114 ! 1115 force(:,:) = force(:,:) - lagrange(i)*dg(:,:,i) 1116 ! 1117 ENDDO 1118 ! 1119 DEALLOCATE( dg_matrix ) 1120 DEALLOCATE( iwork ) 1121 ! 1122 ENDIF 1123 ! 1124#if defined (__DEBUG_CONSTRAINTS) 1125 ! 1126 WRITE( stdout, '(/,5X,"Intermediate forces (Ry/au):",/)') 1127 ! 1128 DO i = 1, nat 1129 ! 1130 WRITE( stdout, '(5X,"atom ",I3," type ",I2,3X,"force = ",3F14.8)' ) & 1131 i, ityp(i), force(:,i) 1132 ! 1133 ENDDO 1134 ! 1135#endif 1136 ! 1137 norm_after = dnrm2( dim, force, 1 ) 1138 ! 1139 IF ( norm_before < norm_after ) THEN 1140 ! 1141 WRITE( stdout, '(/,5X,"norm before = ",F16.10)' ) norm_before 1142 WRITE( stdout, '( 5X,"norm after = ",F16.10)' ) norm_after 1143 ! 1144 CALL errore( 'remove_constr_force', & 1145 'norm(F) before < norm(F) after', 1 ) 1146 ! 1147 ENDIF 1148 ! 1149 DEALLOCATE( dg ) 1150 ! 1151#endif 1152 ! 1153 END SUBROUTINE remove_constr_force 1154 ! 1155 !----------------------------------------------------------------------- 1156 SUBROUTINE remove_constr_vec( nat, tau, & 1157 if_pos, ityp, tau_units, vec ) 1158 !----------------------------------------------------------------------- 1159 ! 1160 ! ... the component of a displacement vector that is orthogonal to the 1161 ! ... ipersurface defined by the constraint equations is removed 1162 ! ... and the corresponding value of the lagrange multiplier computed 1163 ! 1164 IMPLICIT NONE 1165 ! 1166 INTEGER, INTENT(in) :: nat 1167 REAL(DP), INTENT(in) :: tau(:,:) 1168 INTEGER, INTENT(in) :: if_pos(:,:) 1169 INTEGER, INTENT(in) :: ityp(:) 1170 REAL(DP), INTENT(in) :: tau_units 1171 REAL(DP), INTENT(inout) :: vec(:,:) 1172 ! 1173 INTEGER :: i, j, dim 1174 REAL(DP) :: g, ndg, dgidgj 1175 REAL(DP), ALLOCATABLE :: dg(:,:,:), dg_matrix(:,:), lambda(:) 1176 INTEGER, ALLOCATABLE :: iwork(:) 1177 ! 1178 REAL(DP), EXTERNAL :: ddot 1179 ! 1180 ! 1181 dim = 3*nat 1182 ! 1183 ALLOCATE( lambda( nconstr ) ) 1184 ALLOCATE( dg( 3, nat, nconstr ) ) 1185 ! 1186 IF ( nconstr == 1 ) THEN 1187 ! 1188 CALL constraint_grad( 1, nat, tau, & 1189 if_pos, ityp, tau_units, g, dg(:,:,1) ) 1190 ! 1191 lambda(1) = ddot( dim, vec, 1, dg(:,:,1), 1 ) 1192 ! 1193 ndg = ddot( dim, dg(:,:,1), 1, dg(:,:,1), 1 ) 1194 ! 1195 vec(:,:) = vec(:,:) - lambda(1)*dg(:,:,1)/ndg 1196 ! 1197 ELSE 1198 ! 1199 ALLOCATE( dg_matrix( nconstr, nconstr ) ) 1200 ALLOCATE( iwork( nconstr ) ) 1201 ! 1202 DO i = 1, nconstr 1203 ! 1204 CALL constraint_grad( i, nat, tau, & 1205 if_pos, ityp, tau_units, g, dg(:,:,i) ) 1206 ! 1207 ENDDO 1208 ! 1209 DO i = 1, nconstr 1210 ! 1211 dg_matrix(i,i) = ddot( dim, dg(:,:,i), 1, dg(:,:,i), 1 ) 1212 ! 1213 lambda(i) = ddot( dim, vec, 1, dg(:,:,i), 1 ) 1214 ! 1215 DO j = i + 1, nconstr 1216 ! 1217 dgidgj = ddot( dim, dg(:,:,i), 1, dg(:,:,j), 1 ) 1218 ! 1219 dg_matrix(i,j) = dgidgj 1220 dg_matrix(j,i) = dgidgj 1221 ! 1222 ENDDO 1223 ! 1224 ENDDO 1225 ! 1226 CALL DGESV( nconstr, 1, dg_matrix, & 1227 nconstr, iwork, lambda, nconstr, i ) 1228 ! 1229 IF ( i /= 0 ) & 1230 CALL errore( 'remove_constr_vec', & 1231 'error in the solution of the linear system', i ) 1232 ! 1233 DO i = 1, nconstr 1234 ! 1235 vec(:,:) = vec(:,:) - lambda(i)*dg(:,:,i) 1236 ! 1237 ENDDO 1238 ! 1239 DEALLOCATE( dg_matrix ) 1240 DEALLOCATE( iwork ) 1241 ! 1242 ENDIF 1243 ! 1244 DEALLOCATE( lambda, dg ) 1245 ! 1246 END SUBROUTINE remove_constr_vec 1247 ! 1248 !----------------------------------------------------------------------- 1249 SUBROUTINE deallocate_constraint() 1250 !----------------------------------------------------------------------- 1251 ! 1252 IMPLICIT NONE 1253 ! 1254 ! 1255 IF ( allocated( lagrange ) ) DEALLOCATE( lagrange ) 1256 IF ( allocated( constr ) ) DEALLOCATE( constr ) 1257 IF ( allocated( constr_type ) ) DEALLOCATE( constr_type ) 1258 IF ( allocated( constr_target ) ) DEALLOCATE( constr_target ) 1259 IF ( allocated( gp ) ) DEALLOCATE( gp ) 1260 ! 1261 RETURN 1262 ! 1263 END SUBROUTINE deallocate_constraint 1264 ! 1265 !----------------------------------------------------------------------- 1266 FUNCTION cross(A,B) 1267 !----------------------------------------------------------------------- 1268 ! 1269 ! ... cross product 1270 ! 1271 IMPLICIT NONE 1272 ! 1273 REAL(DP),INTENT(in) :: A(3),B(3) 1274 REAL(DP) cross(3) 1275 ! 1276 cross(1) = A(2)*B(3)-A(3)*B(2) 1277 cross(2) = A(3)*B(1)-A(1)*B(3) 1278 cross(3) = A(1)*B(2)-A(2)*B(1) 1279 ! 1280 END FUNCTION 1281 ! 1282 !----------------------------------------------------------------------- 1283 FUNCTION pbc( vect ) 1284 !----------------------------------------------------------------------- 1285 ! 1286 ! ... periodic boundary conditions ( vect is assumed to be given 1287 ! ... in cartesian coordinates and in atomic units ) 1288 ! 1289 USE cell_base, ONLY : at, bg, alat 1290 ! 1291 IMPLICIT NONE 1292 ! 1293 REAL(DP), INTENT(in) :: vect(3) 1294 REAL(DP) :: pbc(3) 1295 ! 1296 ! 1297#if defined (__USE_PBC) 1298 ! 1299 pbc(:) = matmul( vect(:), bg(:,:) )/alat 1300 ! 1301 pbc(:) = pbc(:) - anint( pbc(:) ) 1302 ! 1303 pbc(:) = matmul( at(:,:), pbc(:) )*alat 1304 ! 1305#else 1306 ! 1307 pbc(:) = vect(:) 1308 ! 1309#endif 1310 RETURN 1311 ! 1312 END FUNCTION pbc 1313 ! 1314 !----------------------------------------------------------------------- 1315 SUBROUTINE compute_dmax() 1316 !----------------------------------------------------------------------- 1317 ! 1318 ! ... dmax corresponds to one half the longest diagonal of the cell 1319 ! 1320 USE cell_base, ONLY : at, alat 1321 ! 1322 IMPLICIT NONE 1323 ! 1324 INTEGER :: x,y,z 1325 REAL(DP) :: diago(3) 1326 ! 1327 dmax = 0._dp !norm(at(:,1)+at(:,2)+at(:,3)) 1328 ! 1329 DO z = -1,1,2 1330 DO y = -1,1,2 1331 DO x = -1,1,2 1332 diago = x*at(:,1) + y*at(:,2) + z*at(:,3) 1333 dmax = max(dmax, norm(diago)) 1334 ENDDO 1335 ENDDO 1336 ENDDO 1337 ! 1338 dmax= dmax*alat*.5_dp 1339 ! 1340 RETURN 1341 ! 1342 END SUBROUTINE compute_dmax 1343 ! 1344END MODULE constraints_module 1345