TABLE OF CONTENTS
- ABINIT/glue_lotf
- eval_lotf/eval_Upp_n
- glue_lotf/calc_coord
- glue_lotf/calc_coord_new_d
- glue_lotf/calc_rhop
- glue_lotf/eval_U_n
- glue_lotf/glue_init
- glue_lotf/glue_pair
- glue_lotf/glue_pair_devs
- glue_lotf/rho_devs
- glue_lotf/rhop_value
ABINIT/glue_lotf [ Modules ]
NAME
glue_lotf
FUNCTION
Contains the GLUE procedure and parameters for Lotf
COPYRIGHT
Copyright (C) 2009-2024 ABINIT group (MMancini) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
Parameters for the glue potential ref: Ercolessi et al, Phil Mag A 58(1), 213 (1988)
SOURCE
20 #if defined HAVE_CONFIG_H 21 #include "config.h" 22 #endif 23 24 #include "abi_common.h" 25 26 module glue_lotf 27 28 use defs_basis 29 use m_errors 30 31 implicit none 32 private 33 34 !--Parameters for the glue potential (Ercolessi et al, Phil Mag A 58(1), 213 (1988)) 35 36 !--Function phi_r 37 real(dp),public :: dphi,rcphi 38 real(dp) :: a0I,a1I,a2I,a3I,a4I 39 real(dp) :: a0II,a1II,a2II,a3II,a4II,a5II,a6II 40 41 !--Function rho_r 42 real(dp),public :: rmrho 43 real(dp) :: rbrho 44 real(dp) :: b0I,b1I,b2I,b3I 45 real(dp) :: b0II,b1II,b2II,b3II,b0III,b1III,b2III,b3III 46 47 !--Function U_n 48 real(dp) :: n0U,nsU 49 real(dp) :: c0I,c1I,c2I,c3I,c4I 50 real(dp) :: c0II,c1II,c2II,c3II,c4II,c0III,c1III,c2III,c3III 51 52 !--Auxiliary variables to change glue parameters 53 real(dp),parameter :: ddphi = 0.028779246_dp 54 real(dp),parameter :: scalefactor_phi = one 55 56 57 public :: & 58 glue_init, & 59 glue_pair_devs, & 60 glue_pair, & 61 calc_coord, & 62 calc_rhop, & 63 calc_coord_new_d, & 64 rho_devs, & 65 rhop_value, & 66 eval_U_n, & 67 eval_Upp_n 68 69 contains !===========================================================
eval_lotf/eval_Upp_n [ Functions ]
[ Top ] [ eval_lotf ] [ Functions ]
NAME
eval_Upp_n
FUNCTION
INPUTS
SOURCE
582 subroutine eval_Upp_n(coordatom_i,up_dum,upp_dum) 583 584 implicit none 585 586 !Arguments ------------------------ 587 real(dp),intent(in) :: coordatom_i 588 real(dp),intent(out) :: up_dum,upp_dum 589 !Local --------------------------- 590 real(dp) :: cns,cns2,cns3,cn0,cn02,cn03 591 character(len=500) :: msg 592 593 ! ************************************************************************* 594 595 !--Evaluate U and Up for atom i 596 up_dum = zero 597 upp_dum = zero 598 599 if (coordatom_i < zero) then 600 601 write(msg,'(3a,i8,2a)')& 602 & 'ERROR: negative coordination!!! ',ch10,& 603 & 'coord = ', coordatom_i,ch10,& 604 & 'The coordination cannot be negative!' 605 ABI_ERROR(msg) 606 607 elseif (coordatom_i < nsU) then 608 609 cns = coordatom_i - nsU 610 cns2 = cns * cns 611 cns3 = cns2 * cns 612 613 up_dum = 4.d0*c4I*cns3 + 3.d0*c3I*cns2 + 2.d0*c2I*cns + c1I 614 upp_dum = 12.d0*c4I*cns2 + 6.d0*c3I*cns + 2.d0*c2I 615 616 elseif (coordatom_i < n0U) then 617 618 cn0 = coordatom_i - n0U 619 cn02 = cn0 * cn0 620 cn03 = cn02 * cn0 621 622 up_dum = 4.d0*c4II*cn03 + 3.d0*c3II*cn02 + 2.d0*c2II*cn0 + c1II 623 upp_dum = 12.d0*c4II*cn02 + 6.d0*c3II*cn0 + 2.d0*c2II 624 625 elseif (coordatom_i >= n0U) then 626 627 cn0 = coordatom_i - n0U 628 cn02 = cn0 * cn0 629 630 up_dum = 3.d0*c3III*cn02 + 2.d0*c2III*cn0 + c1III 631 upp_dum = 6.d0*c3III*cn0 + 2.d0*c2III 632 633 endif 634 end subroutine eval_Upp_n 635 636 end module glue_lotf
glue_lotf/calc_coord [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
calc_coord
FUNCTION
INPUTS
SOURCE
304 subroutine calc_coord(r_au,coordatom_dum) 305 306 implicit none 307 308 !Arguments ------------------------ 309 real(dp) :: r_au,coordatom_dum 310 !Local --------------------------- 311 real(dp) :: rst,r,rho_neigh 312 ! ********************************************************************* 313 ! --Adds the density of a single given neighbour 314 315 rst = sqrt(r_au) 316 rho_neigh = zero 317 318 if (rst < dphi) then 319 r = rst - dphi 320 rho_neigh = b3I*r*r*r + b2I*r*r + b1I*r + b0I 321 elseif (rst < rbrho) then 322 r = rst - dphi 323 rho_neigh = b3II*r*r*r + b2II*r*r + b1II*r + b0II 324 elseif (rst < rmrho) then 325 r = rst - rmrho 326 rho_neigh = b3III*r*r*r + b2III*r*r + b1III*r + b0III 327 end if 328 329 coordatom_dum = coordatom_dum + rho_neigh 330 331 end subroutine calc_coord
glue_lotf/calc_coord_new_d [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
calc_coord_new_d
FUNCTION
INPUTS
SOURCE
386 subroutine calc_coord_new_d(r_au,alpha_d,coordatom_dum) 387 388 implicit none 389 390 !Arguments ------------------------ 391 real(dp),intent(in) :: r_au,alpha_d 392 real(dp),intent(out) ::coordatom_dum 393 !Local --------------------------- 394 real(dp) :: rst,r,rho_neigh 395 ! ********************************************************************* 396 397 ! --Adds the density of a single given neighbour 398 rst = sqrt(r_au) 399 400 rho_neigh = zero 401 402 if (rst < alpha_d) then 403 r = rst - alpha_d 404 rho_neigh = b3I*r*r*r + b2I*r*r + b1I*r + b0I 405 elseif (rst < rbrho) then 406 r = rst - alpha_d 407 rho_neigh = b3II*r*r*r + b2II*r*r + b1II*r + b0II 408 409 ! --Restriction on rmrho to keep rho and rhop continuous 410 elseif (rst < (rmrho+alpha_d-dphi)) then 411 r = rst - (rmrho+alpha_d-dphi) 412 rho_neigh = b3III*r*r*r + b2III*r*r + b1III*r + b0III 413 end if 414 415 coordatom_dum = coordatom_dum + rho_neigh 416 end subroutine calc_coord_new_d
glue_lotf/calc_rhop [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
calc_rhop
FUNCTION
INPUTS
SOURCE
346 subroutine calc_rhop(r_st,rhop_dum) 347 348 implicit none 349 350 ! Arguments ------------------------ 351 real(dp),intent(in) :: r_st 352 real(dp),intent(out) :: rhop_dum 353 ! Local --------------------------- 354 real(dp) :: r,r2 355 ! ********************************************************************* 356 357 rhop_dum = zero 358 359 if (r_st < dphi) then 360 r = r_st - dphi 361 r2 = r*r 362 rhop_dum = (3.d0*b3I*r2 + 2.d0*b2I*r + b1I) 363 elseif (r_st < rbrho) then 364 r = r_st - dphi 365 r2 = r*r 366 rhop_dum = (3.d0*b3II*r2 + 2.d0*b2II*r + b1II) 367 elseif (r_st < rmrho) then 368 r = r_st - rmrho 369 r2 = r*r 370 rhop_dum = (3.d0*b3III*r2 + 2.d0*b2III*r + b1III) 371 end if 372 end subroutine calc_rhop
glue_lotf/eval_U_n [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
eval_U_n
FUNCTION
INPUTS
SOURCE
516 subroutine eval_U_n(coordatom_i,epot_dum2,up_dum) 517 518 implicit none 519 520 !Arguments ------------------------ 521 real(dp),intent(in) :: coordatom_i 522 real(dp),intent(out) :: epot_dum2,up_dum 523 !Local --------------------------- 524 real(dp) :: cns,cns2,cns3,cns4,cn0,cn02,cn03,cn04 525 character(len=500) :: msg 526 527 !--Evaluate U and Up for atom i 528 epot_dum2 = zero 529 up_dum = zero 530 531 if (coordatom_i < zero) then 532 533 write(msg,'(3a,i8,2a)')& 534 & 'ERROR: negative coordination!!! ',ch10,& 535 & 'coord = ', coordatom_i,ch10,& 536 & 'The coordination cannot be negative!' 537 ABI_ERROR(msg) 538 539 elseif (coordatom_i < nsU) then 540 541 cns = coordatom_i - nsU 542 cns2 = cns * cns 543 cns3 = cns2 * cns 544 cns4 = cns3 * cns 545 546 epot_dum2 = c4I*cns4 + c3I*cns3 + c2I*cns2 + c1I*cns + c0I 547 up_dum = 4.d0*c4I*cns3 + 3.d0*c3I*cns2 + 2.d0*c2I*cns + c1I 548 549 elseif (coordatom_i < n0U) then 550 551 cn0 = coordatom_i - n0U 552 cn02 = cn0 * cn0 553 cn03 = cn02 * cn0 554 cn04 = cn03 * cn0 555 556 epot_dum2 = c4II*cn04 + c3II*cn03 + c2II*cn02 + c1II*cn0 + c0II 557 up_dum = 4.d0*c4II*cn03 + 3.d0*c3II*cn02 + 2.d0*c2II*cn0 + c1II 558 559 elseif (coordatom_i >= n0U) then 560 561 cn0 = coordatom_i - n0U 562 cn02 = cn0 * cn0 563 cn03 = cn02 * cn0 564 565 epot_dum2 = c3III*cn03 + c2III*cn02 + c1III*cn0 + c0III 566 up_dum = 3.d0*c3III*cn02 + 2.d0*c2III*cn0 + c1III 567 568 endif 569 end subroutine eval_U_n
glue_lotf/glue_init [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
glue_init
FUNCTION
INPUTS
SOURCE
81 subroutine glue_init() 82 83 !--Function phi_r 84 ! ************************************************************************* 85 dphi = (0.2878207442141723d+01 + ddphi) / 0.529177d0 86 rcphi = (0.3700000000000000d+01 + ddphi) / 0.529177d0 87 88 a0I = -0.8000000000000000d-01 / (13.6057d0 * 2.d0) *& 89 scalefactor_phi 90 a1I = 0.0000000000000000d+00 / (13.6057d0 * 2.d0) * (0.529177d0) *& 91 scalefactor_phi 92 a2I = 0.7619231375231362d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**2 *& 93 scalefactor_phi 94 a3I = -0.8333333333333333d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**3 *& 95 scalefactor_phi 96 a4I = -0.1211483464993159d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**4 *& 97 scalefactor_phi 98 99 ! a0II = -0.8000000000000000d-01 / (13.6057d0 * 2.d0) 100 a0II = a0I 101 ! a1II = 0.0000000000000000d+00 / (13.6057d0 * 2.d0) * (0.529177d0) 102 a1II = a1I 103 ! a2II = 0.7619231375231362d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**2 104 a2II = a2I 105 a3II = -0.8333333333333333d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**3 *& 106 scalefactor_phi 107 a4II = -0.1096009851140349d+01 / (13.6057d0 * 2.d0) * (0.529177d0)**4 *& 108 scalefactor_phi 109 a5II = 0.2158417178555998d+01 / (13.6057d0 * 2.d0) * (0.529177d0)**5 *& 110 scalefactor_phi 111 a6II = -0.9128915709636862d+00 / (13.6057d0 * 2.d0) * (0.529177d0)**6 *& 112 scalefactor_phi 113 114 !--Function rho_r 115 rbrho = 0.3500000000000000d+01 / 0.529177d0 116 rmrho = (0.3900000000000000d+01 + ddphi) / 0.529177d0 117 118 b0I = 0.1000000000000000d+01 119 b1I = -0.6800000000000000d+00 * (0.529177d0) 120 b2I = 0.7500000000000000d+00 * (0.529177d0)**2 121 b3I = -0.1333333333333333d+01 * (0.529177d0)**3 122 123 ! b0II = 0.1000000000000000d+01 124 b0II = b0I 125 ! b1II = -0.6800000000000000d+00 * (0.529177d0) 126 b1II = b1I 127 ! b2II = 0.7500000000000000d+00 * (0.529177d0)**2 128 b2II = b2I 129 b3II = -0.1527241171296038d+01 * (0.529177d0)**3 130 131 b0III = 0.0000000000000000d+00 132 b1III = 0.0000000000000000d+00 * (0.529177d0) 133 b2III = 0.5578188675490974d+01 * (0.529177d0)**2 134 b3III = 0.6132971688727435d+01 * (0.529177d0)**3 135 136 !--Function U_n 137 n0U = 0.1200000000000000d+02 138 nsU = 0.9358157767784574d+01 139 140 c0I = -0.2793388616771698d+01 / (13.6057d0 * 2.d0) * scalefactor_phi 141 c1I = -0.3419999999999999d+00 / (13.6057d0 * 2.d0) * scalefactor_phi 142 c2I = 0.3902327808424106d-01 / (13.6057d0 * 2.d0) * scalefactor_phi 143 c3I = 0.7558829951858879d-02 / (13.6057d0 * 2.d0) * scalefactor_phi 144 c4I = 0.3090472511796849d-03 / (13.6057d0 * 2.d0) * scalefactor_phi 145 146 c0II = -0.3300000000000000d+01 / (13.6057d0 * 2.d0) * scalefactor_phi 147 c1II = 0.0000000000000000d+00 / (13.6057d0 * 2.d0) * scalefactor_phi 148 c2II = 0.8618226772941980d-01 / (13.6057d0 * 2.d0) * scalefactor_phi 149 c3II = 0.4341701445034724d-02 / (13.6057d0 * 2.d0) * scalefactor_phi 150 c4II = -0.3044398779375916d-03 / (13.6057d0 * 2.d0) * scalefactor_phi 151 152 ! c0III = -0.3300000000000000d+01 / (13.6057d0 * 2.d0) 153 c0III = c0II 154 ! c1III = 0.0000000000000000d+00 / (13.6057d0 * 2.d0) 155 c1III = c1II 156 ! c2III = 0.8618226772941980d-01 / (13.6057d0 * 2.d0) 157 c2III = c2II 158 c3III = 0.4325981467602070d-02 / (13.6057d0 * 2.d0) * scalefactor_phi 159 160 end subroutine glue_init
glue_lotf/glue_pair [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
glue_pair
FUNCTION
INPUTS
SOURCE
248 subroutine glue_pair(RD,r_au,epot_2,fdum) 249 250 implicit none 251 252 !Arguments ------------------------ 253 real(dp) :: RD(3), fdum(3), r_au, epot_2 254 !Local --------------------------- 255 real(dp) :: rst, r, r2, r3, r4, r5, r6 256 integer :: k 257 ! ********************************************************************* 258 ! --The calculation for a given pair 259 rst = sqrt(r_au) 260 r = rst - dphi 261 r2 = r*r 262 r3 = r2*r 263 r4 = r3*r 264 265 ! --Energy 266 epot_2 = zero 267 if (rst < dphi) then 268 epot_2 = a4I*r4 + a3I*r3 + a2I*r2 + a1I*r + a0I 269 elseif (rst < rcphi) then 270 r5 = r4*r 271 r6 = r5*r 272 epot_2 = a6II*r6 + a5II*r5 + a4II*r4 + a3II*r3 + a2II*r2 +& 273 a1II*r + a0II 274 end if 275 276 ! --Forces (without the minus sign, because we use rji) 277 fdum(:) = zero 278 279 do k =1,3 280 if (rst < dphi) then 281 fdum(k) = ( 4.d0*a4I*r3 + 3.d0*a3I*r2 + 2.d0*a2I*r + a1I ) * & 282 RD(k)/rst 283 elseif (rst < rcphi) then 284 fdum(k) = ( 6.d0*a6II*r5 + 5.d0*a5II*r4 + 4.d0*a4II*r3 + & 285 3.d0*a3II*r2 + 2.d0*a2II*r + a1II ) * & 286 RD(k)/rst 287 end if 288 end do 289 290 end subroutine glue_pair
glue_lotf/glue_pair_devs [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
glue_pair_devs
FUNCTION
INPUTS
SOURCE
173 subroutine glue_pair_devs(alpha_dum,RD,r_au,epot_2,fdum,dfdum) 174 175 implicit none 176 177 !Arguments ------------------------ 178 real(dp),intent(in) :: alpha_dum(3),RD(3) 179 real(dp),intent(in) :: r_au 180 real(dp),intent(out) :: epot_2 181 real(dp),intent(out) :: fdum(3),dfdum(3,2) 182 !Local --------------------------- 183 real(dp) :: rst, r, r2, r3, r4, r5, r6 184 real(dp) :: drcphi, dphi2, rcphi2, scale_factor 185 ! ********************************************************************* 186 187 ! --The calculation for a given pair 188 drcphi = rcphi - dphi 189 dphi2 = alpha_dum(1) 190 191 ! --Rewritten 192 rcphi2 = dphi2 + drcphi 193 194 scale_factor = alpha_dum(2) 195 196 rst = sqrt(r_au) 197 198 r = rst - dphi2 199 r2 = r*r 200 r3 = r2*r 201 r4 = r3*r 202 203 ! --Energy 204 epot_2 = zero 205 206 if (rst < dphi2) then 207 epot_2 = scale_factor * (a4I*r4 + a3I*r3 + a2I*r2 + a1I*r + a0I) 208 elseif (rst < rcphi2) then 209 r5 = r4*r 210 r6 = r5*r 211 epot_2 = scale_factor * (a6II*r6 + a5II*r5 + a4II*r4 + a3II*r3 + a2II*r2 +& 212 a1II*r + a0II) 213 end if 214 215 ! --Forces (without the minus sign, because we use rji) and derivatives w.r.t. parameter d 216 fdum(:) = zero 217 dfdum(:,:) = zero 218 219 if (rst < dphi2) then 220 fdum(:) = scale_factor * (4.d0*a4I*r3 + 3.d0*a3I*r2 + 2.d0*a2I*r + a1I) * RD(:)/rst 221 dfdum(:,1) = -scale_factor * (12.d0*a4I*r2 + 6.d0*a3I*r + 2.d0*a2I) * RD(:)/rst 222 dfdum(:,2) = fdum(:) / scale_factor 223 224 elseif (rst < rcphi2) then 225 fdum(:) = (scale_factor) * & 226 & (6.d0*a6II*r5 + 5.d0*a5II*r4 + 4.d0*a4II*r3 + 3.d0*a3II*r2 + 2.d0*a2II*r + a1II) * & 227 & RD(:)/rst 228 229 dfdum(:,1) = -scale_factor * & 230 & (30.d0*a6II*r4 + 20.d0*a5II*r3 + 12.d0*a4II*r2 + 6.d0*a3II*r + 2.d0*a2II) * RD(:)/rst 231 232 dfdum(:,2) = fdum(:) / scale_factor 233 234 end if 235 end subroutine glue_pair_devs
glue_lotf/rho_devs [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
rho_devs
FUNCTION
INPUTS
SOURCE
429 subroutine rho_devs(r_au,alpha_d,rho_neigh_d,rho_neigh_p,rho_neigh_pd) 430 431 real(dp),intent(in) :: r_au,alpha_d 432 real(dp),intent(out) :: rho_neigh_d,rho_neigh_p,rho_neigh_pd 433 ! Local --------------------------- 434 real(dp) :: rst,r 435 ! ********************************************************************* 436 437 ! --Adds the density of a single given neighbour 438 rst = sqrt(r_au) 439 440 rho_neigh_p = zero 441 rho_neigh_d = zero 442 rho_neigh_pd = zero 443 444 if (rst < alpha_d) then 445 r = rst - alpha_d 446 rho_neigh_d = - ( 3.d0*b3I*r*r + 2.d0*b2I*r + b1I ) 447 rho_neigh_p = - rho_neigh_d 448 rho_neigh_pd = - ( 6.d0*b3I*r + 2.d0*b2I ) 449 450 elseif (rst < rbrho) then 451 r = rst - alpha_d 452 rho_neigh_d = - ( 3.d0*b3II*r*r + 2.d0*b2II*r + b1II ) 453 rho_neigh_p = - rho_neigh_d 454 rho_neigh_pd = - ( 6.d0*b3II*r + 2.d0*b2II ) 455 456 ! --Restriction on rmrho to keep rho and rhop continuous 457 elseif (rst < (rmrho+alpha_d-dphi)) then 458 r = rst - (rmrho+alpha_d-dphi) 459 rho_neigh_d = zero 460 rho_neigh_p = 3.d0*b3III*r*r + 2.d0*b2III*r + b1III 461 rho_neigh_pd = zero 462 463 end if 464 end subroutine rho_devs
glue_lotf/rhop_value [ Functions ]
[ Top ] [ glue_lotf ] [ Functions ]
NAME
rhop_value
FUNCTION
INPUTS
SOURCE
478 subroutine rhop_value(rst,alpha_d,rhop_dum) 479 480 implicit none 481 482 !Arguments ------------------------ 483 real(dp),intent(in) :: rst,alpha_d 484 real(dp),intent(out) :: rhop_dum 485 !Local --------------------------- 486 real(dp) :: r 487 ! ********************************************************************* 488 489 ! --Adds the density of a single given neighbour 490 rhop_dum = zero 491 492 if (rst < alpha_d) then 493 r = rst - alpha_d 494 rhop_dum = 3.d0*b3I*r*r + 2.d0*b2I*r + b1I 495 elseif (rst < rbrho) then 496 r = rst - alpha_d 497 rhop_dum = 3.d0*b3II*r*r + 2.d0*b2II*r + b1II 498 ! --Restriction on rmrho to keep rho and rhop continuous 499 elseif (rst < (rmrho+alpha_d-dphi)) then 500 r = rst - (rmrho+alpha_d-dphi) 501 rhop_dum = 3.d0*b3III*r*r + 2.d0*b2III*r + b1III 502 end if 503 end subroutine rhop_value