1! 2! CDDL HEADER START 3! 4! The contents of this file are subject to the terms of the Common Development 5! and Distribution License Version 1.0 (the "License"). 6! 7! You can obtain a copy of the license at 8! http://www.opensource.org/licenses/CDDL-1.0. See the License for the 9! specific language governing permissions and limitations under the License. 10! 11! When distributing Covered Code, include this CDDL HEADER in each file and 12! include the License file in a prominent location with the name LICENSE.CDDL. 13! If applicable, add the following below this CDDL HEADER, with the fields 14! enclosed by brackets "[]" replaced with your own identifying information: 15! 16! Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved. 17! 18! CDDL HEADER END 19! 20 21! 22! Copyright (c) 2013, Regents of the University of Minnesota. 23! All rights reserved. 24! 25! Contributors: 26! Ryan S. Elliott 27! Ellad B. Tadmor 28! Valeriu Smirichinski 29! Stephen M. Whalen 30! 31 32!**************************************************************************** 33!** 34!** MODULE EAM_ErcolessiAdams_1994_Al 35!** 36!** Ercolessi-Adams pair functional model for Al 37!** 38!** Reference: F. Ercolessi and J. B. Adams, Europhys. Lett. 26, 583 (1994). 39!** http://www.sissa.it/furio/potentials/Al/ 40!** 41!** Language: Fortran 2003 42!** 43!** 44!**************************************************************************** 45module EAM_ErcolessiAdams_1994_Al 46 47use, intrinsic :: iso_c_binding 48use kim_model_headers_module 49implicit none 50 51save 52private 53public Compute_Energy_Forces, & 54 destroy, & 55 compute_arguments_create, & 56 compute_arguments_destroy, & 57 model_cutoff, & 58 speccode, & 59 buffer_type 60 61! Below are the definitions and values of all Model parameters 62integer(c_int), parameter :: cd = c_double ! used for literal constants 63integer(c_int), parameter :: DIM = 3 ! dimensionality of space 64integer(c_int), parameter :: speccode = 1 ! internal species code 65real(c_double), parameter :: model_cutoff = 5.55805441821810_cd ! cutoff radius 66 ! in angstroms 67real(c_double), parameter :: model_cutsq = model_cutoff**2 68 69!------------------------------------------------------------------------------- 70! Below are the definitions and values of all additional model parameters 71! 72! Recall that the Fortran 2003 format for declaring parameters is as follows: 73! 74integer(c_int), parameter :: npt= 17 75integer(c_int), parameter :: nuu= 13 76real(c_double), parameter :: rhomin = 0.0_cd 77 78! pair potential data 79! 80real(c_double), parameter :: x(npt) = (/ .202111069753385E+01_cd, & 81 .227374953472558E+01_cd, & 82 .252638837191732E+01_cd, & 83 .277902720910905E+01_cd, & 84 .303166604630078E+01_cd, & 85 .328430488349251E+01_cd, & 86 .353694372068424E+01_cd, & 87 .378958255787597E+01_cd, & 88 .404222139506771E+01_cd, & 89 .429486023225944E+01_cd, & 90 .454749906945117E+01_cd, & 91 .480013790664290E+01_cd, & 92 .505277674383463E+01_cd, & 93 .530541558102636E+01_cd, & 94 .555805441821810E+01_cd, & 95 .555807968210182E+01_cd, & 96 .555810494598553E+01_cd /) 97 98real(c_double), parameter :: yv2(npt) = (/ .196016472197158E+01_cd, & 99 .682724240745344E+00_cd, & 100 .147370824539188E+00_cd, & 101 -.188188235860390E-01_cd, & 102 -.576011902692490E-01_cd, & 103 -.519846499644276E-01_cd, & 104 -.376352484845919E-01_cd, & 105 -.373737879689433E-01_cd, & 106 -.531351030124350E-01_cd, & 107 -.632864983555742E-01_cd, & 108 -.548103623840369E-01_cd, & 109 -.372889232343935E-01_cd, & 110 -.188876517630154E-01_cd, & 111 -.585239362533525E-02_cd, & 112 .000000000000000E+00_cd, & 113 .000000000000000E+00_cd, & 114 .000000000000000E+00_cd /) 115 116real(c_double), parameter :: Bv2(npt) = (/ -.702739315585347E+01_cd, & 117 -.333140549270729E+01_cd, & 118 -.117329394261502E+01_cd, & 119 -.306003283486901E+00_cd, & 120 -.366656699104026E-01_cd, & 121 .588330899204400E-01_cd, & 122 .384220572312032E-01_cd, & 123 -.390223173707191E-01_cd, & 124 -.663882722510521E-01_cd, & 125 -.312918894386669E-02_cd, & 126 .590118945294245E-01_cd, & 127 .757939459148246E-01_cd, & 128 .643822548468606E-01_cd, & 129 .399750987463792E-01_cd, & 130 .177103852679117E-05_cd, & 131 -.590423369301474E-06_cd, & 132 .590654950414731E-06_cd /) 133 134real(c_double), parameter :: Cv2(npt) = (/ .877545959718548E+01_cd, & 135 .585407125495837E+01_cd, & 136 .268820820643116E+01_cd, & 137 .744718689404422E+00_cd, & 138 .321378734769888E+00_cd, & 139 .566263292669091E-01_cd, & 140 -.137417679148505E+00_cd, & 141 -.169124163201523E+00_cd, & 142 .608037039066423E-01_cd, & 143 .189589640245655E+00_cd, & 144 .563784150384640E-01_cd, & 145 .100486298765028E-01_cd, & 146 -.552186092621482E-01_cd, & 147 -.413902746758285E-01_cd, & 148 -.116832934994489E+00_cd, & 149 .233610871054729E-01_cd, & 150 .233885865725971E-01_cd /) 151 152real(c_double), parameter :: Dv2(npt) = (/ -.385449887634130E+01_cd, & 153 -.417706040200591E+01_cd, & 154 -.256425277368288E+01_cd, & 155 -.558557503589276E+00_cd, & 156 -.349316054551627E+00_cd, & 157 -.256022933201611E+00_cd, & 158 -.418337423301704E-01_cd, & 159 .303368330939646E+00_cd, & 160 .169921006301015E+00_cd, & 161 -.175759761362548E+00_cd, & 162 -.611278214082881E-01_cd, & 163 -.861140219824535E-01_cd, & 164 .182451950513387E-01_cd, & 165 -.995395392057973E-01_cd, & 166 .184972909229936E+04_cd, & 167 .362829766922787E+00_cd, & 168 .362829766922787E+00_cd /) 169 170! electron density data 171! 172real(c_double), parameter :: yrh(npt) = (/ .865674623712589E-01_cd, & 173 .925214702944478E-01_cd, & 174 .862003123832002E-01_cd, & 175 .762736292751052E-01_cd, & 176 .606481841271735E-01_cd, & 177 .466030959588197E-01_cd, & 178 .338740138848363E-01_cd, & 179 .232572661705343E-01_cd, & 180 .109046405489829E-01_cd, & 181 .524910605677597E-02_cd, & 182 .391702419142291E-02_cd, & 183 .308277776293383E-02_cd, & 184 .250214745349505E-02_cd, & 185 .147220513798186E-02_cd, & 186 .000000000000000E+00_cd, & 187 .000000000000000E+00_cd, & 188 .000000000000000E+00_cd /) 189 190real(c_double), parameter :: Brh(npt) = (/ .608555214104682E-01_cd, & 191 -.800158928716306E-02_cd, & 192 -.332089451111092E-01_cd, & 193 -.521001991705069E-01_cd, & 194 -.618130637429111E-01_cd, & 195 -.529750064268036E-01_cd, & 196 -.442210477548108E-01_cd, & 197 -.473645664984640E-01_cd, & 198 -.390741582571631E-01_cd, & 199 -.101795580610560E-01_cd, & 200 -.318316981110289E-02_cd, & 201 -.281217210746153E-02_cd, & 202 -.236932031483360E-02_cd, & 203 -.683554708271547E-02_cd, & 204 -.638718204858808E-06_cd, & 205 .212925486831149E-06_cd, & 206 -.212983742465787E-06_cd /) 207 208real(c_double), parameter :: Crh(npt) = (/ -.170233687052940E+00_cd, & 209 -.102317878901959E+00_cd, & 210 .254162872544396E-02_cd, & 211 -.773173610292656E-01_cd, & 212 .388717099948882E-01_cd, & 213 -.388873819867093E-02_cd, & 214 .385388290924526E-01_cd, & 215 -.509815666327127E-01_cd, & 216 .837968231208082E-01_cd, & 217 .305743500420042E-01_cd, & 218 -.288110886134041E-02_cd, & 219 .434959924771674E-02_cd, & 220 -.259669459714693E-02_cd, & 221 -.150816117849093E-01_cd, & 222 .421356801161513E-01_cd, & 223 -.842575249165724E-02_cd, & 224 -.843267014952237E-02_cd /) 225 226real(c_double), parameter :: Drh(npt) = (/ .896085612514625E-01_cd, & 227 .138352319847830E+00_cd, & 228 -.105366473134009E+00_cd, & 229 .153300619856764E+00_cd, & 230 -.564184148788224E-01_cd, & 231 .559792096400504E-01_cd, & 232 -.118113795329664E+00_cd, & 233 .177827488509794E+00_cd, & 234 -.702220789044304E-01_cd, & 235 -.441413511810337E-01_cd, & 236 .954024354744484E-02_cd, & 237 -.916498550800407E-02_cd, & 238 -.164726813535368E-01_cd, & 239 .754928689733184E-01_cd, & 240 -.667110847110954E+03_cd, & 241 -.912720300911022E-01_cd, & 242 -.912720300911022E-01_cd /) 243 244! Embedding function data 245! 246real(c_double), parameter :: xuu(nuu) = (/ .000000000000000E+00_cd, & 247 .100000000000000E+00_cd, & 248 .200000000000000E+00_cd, & 249 .300000000000000E+00_cd, & 250 .400000000000000E+00_cd, & 251 .500000000000000E+00_cd, & 252 .600000000000000E+00_cd, & 253 .700000000000000E+00_cd, & 254 .800000000000000E+00_cd, & 255 .900000000000000E+00_cd, & 256 .100000000000000E+01_cd, & 257 .110000000000000E+01_cd, & 258 .120000000000000E+01_cd /) 259 260real(c_double), parameter :: yuu(nuu) = (/ .000000000000000E+00_cd, & 261 -.113953324143752E+01_cd, & 262 -.145709859805864E+01_cd, & 263 -.174913308002738E+01_cd, & 264 -.202960322136630E+01_cd, & 265 -.225202324967546E+01_cd, & 266 -.242723053979436E+01_cd, & 267 -.255171976467357E+01_cd, & 268 -.260521638832322E+01_cd, & 269 -.264397894381693E+01_cd, & 270 -.265707884842034E+01_cd, & 271 -.264564149400021E+01_cd, & 272 -.260870604452106E+01_cd /) 273 274real(c_double), parameter :: Buu(nuu) = (/ -.183757286015853E+02_cd, & 275 -.574233124410516E+01_cd, & 276 -.236790436375322E+01_cd, & 277 -.307404645857774E+01_cd, & 278 -.251104850116555E+01_cd, & 279 -.196846462620234E+01_cd, & 280 -.154391254686695E+01_cd, & 281 -.846780636273251E+00_cd, & 282 -.408540363905760E+00_cd, & 283 -.286833282404628E+00_cd, & 284 -.309389414590161E-06_cd, & 285 .236958014464143E+00_cd, & 286 .503352368511243E+00_cd /) 287 288real(c_double), parameter :: Cuu(nuu) = (/ .830779120415016E+02_cd, & 289 .432560615333001E+02_cd, & 290 -.951179272978074E+01_cd, & 291 .245037178153561E+01_cd, & 292 .317960779258630E+01_cd, & 293 .224623095704576E+01_cd, & 294 .199928983630817E+01_cd, & 295 .497202926962879E+01_cd, & 296 -.589626545953876E+00_cd, & 297 .180669736096520E+01_cd, & 298 .106163236918694E+01_cd, & 299 .130795086934864E+01_cd, & 300 .135599267112235E+01_cd /) 301 302real(c_double), parameter :: Duu(nuu) = (/ -.132739501694005E+03_cd, & 303 -.175892847543603E+03_cd, & 304 .398738817043878E+02_cd, & 305 .243078670350231E+01_cd, & 306 -.311125611846847E+01_cd, & 307 -.823137069125319E+00_cd, & 308 .990913144440207E+01_cd, & 309 -.185388527186089E+02_cd, & 310 .798774635639692E+01_cd, & 311 -.248354997259420E+01_cd, & 312 .821061667205675E+00_cd, & 313 .160139339245701E+00_cd, & 314 .160139339245701E+00_cd /) 315 316! Buffer 317! 318! (irlast is the last entry into the pair potential and electron density 319! spline tables. It is stored to speed calculations since the next value 320! of r is likely close to the last one.) 321! 322! (ielast is the last entry into the embedding function spline table. 323! It is stored to speed calculations since the next value of rho is 324! likely close to the last one.) 325 326type, bind(c) :: buffer_type 327 real(c_double) :: influence_distance 328 real(c_double) :: cutoff(1) 329 integer(c_int) :: irlast 330 integer(c_int) :: ielast 331 integer(c_int) :: & 332 model_will_not_request_neighbors_of_noncontributing_particles(1) 333end type buffer_type 334 335contains 336 337!------------------------------------------------------------------------------- 338! 339! Calculate pair potential phi(r) 340! 341!------------------------------------------------------------------------------- 342! The "recursive" keyword is added below make this routine thread safe 343recursive subroutine calc_phi(r,phi,irlast) 344 implicit none 345 346 !-- Transferred variables 347 real(c_double), intent(in) :: r 348 real(c_double), intent(out) :: phi 349 integer(c_int), intent(inout) :: irlast 350 351 !-- Local variables 352 integer(c_int) i 353 real(c_double) dx 354 355 if (r .gt. model_cutoff) then 356 ! Argument exceeds cutoff radius 357 phi = 0.0_cd 358 else 359 i = seval_i(npt,r,x,irlast) 360 dx = r - x(i) 361 phi = yv2(i) + dx*(Bv2(i) + dx*(Cv2(i) + dx*Dv2(i))) 362 endif 363 364 return 365end subroutine calc_phi 366 367!------------------------------------------------------------------------------- 368! 369! Calculate pair potential phi(r) and its derivative dphi(r) 370! 371!------------------------------------------------------------------------------- 372! The "recursive" keyword is added below make this routine thread safe 373recursive subroutine calc_phi_dphi(r,phi,dphi,irlast) 374 implicit none 375 376 !-- Transferred variables 377 real(c_double), intent(in) :: r 378 real(c_double), intent(out) :: phi,dphi 379 integer(c_int), intent(inout) :: irlast 380 381 !-- Local variables 382 integer(c_int) i 383 real(c_double) dx 384 385 if (r .gt. model_cutoff) then 386 ! Argument exceeds cutoff radius 387 phi = 0.0_cd 388 dphi = 0.0_cd 389 else 390 i = seval_i(npt,r,x,irlast) 391 dx = r - x(i) 392 phi = yv2(i) + dx*(Bv2(i) + dx*(Cv2(i) + dx*Dv2(i))) 393 dphi = Bv2(i) + dx*(2.0_cd*Cv2(i) + 3.0_cd*dx*Dv2(i)) 394 endif 395 396 return 397end subroutine calc_phi_dphi 398 399!------------------------------------------------------------------------------- 400! 401! Calculate electron density g(r) 402! 403!------------------------------------------------------------------------------- 404! The "recursive" keyword is added below make this routine thread safe 405recursive subroutine calc_g(r,g,irlast) 406 implicit none 407 408 !-- Transferred variables 409 real(c_double), intent(in) :: r 410 real(c_double), intent(out) :: g 411 integer(c_int), intent(inout) :: irlast 412 413 !-- Local variables 414 integer(c_int) i 415 real(c_double) dx 416 417 if (r .gt. model_cutoff) then 418 ! Argument exceeds cutoff radius 419 g = 0.0_cd 420 else 421 i = seval_i(npt,r,x,irlast) 422 dx = r - x(i) 423 g = yrh(i) + dx*(Brh(i) + dx*(Crh(i) + dx*Drh(i))) 424 endif 425 426 return 427end subroutine calc_g 428 429!------------------------------------------------------------------------------- 430! 431! Calculate electron density derivative dg(r) 432! 433!------------------------------------------------------------------------------- 434! The "recursive" keyword is added below make this routine thread safe 435recursive subroutine calc_dg(r,dg,irlast) 436 implicit none 437 438 !-- Transferred variables 439 real(c_double), intent(in) :: r 440 real(c_double), intent(out) :: dg 441 integer(c_int), intent(inout) :: irlast 442 443 !-- Local variables 444 integer(c_int) i 445 real(c_double) dx 446 447 if (r .gt. model_cutoff) then 448 ! Argument exceeds cutoff radius 449 dg = 0.0_cd 450 else 451 i = seval_i(npt,r,x,irlast) 452 dx = r - x(i) 453 dg = Brh(i) + dx*(2.0_cd*Crh(i) + 3.0_cd*dx*Drh(i)) 454 endif 455 456 return 457end subroutine calc_dg 458 459!------------------------------------------------------------------------------- 460! 461! Calculate embedding function U(rho) 462! 463!------------------------------------------------------------------------------- 464! The "recursive" keyword is added below make this routine thread safe 465recursive subroutine calc_U(rho,U,ielast) 466 implicit none 467 468 !-- Transferred variables 469 real(c_double), intent(in) :: rho 470 real(c_double), intent(out) :: U 471 integer(c_int), intent(inout) :: ielast 472 473 !-- Local variables 474 integer(c_int) i 475 real(c_double) dx 476 477 if (rho .le. rhomin) then 478 ! Argument less than the minimum stored value 479 U = 0.0_cd 480 else 481 i = seval_i(nuu,rho,xuu,ielast) 482 dx = rho - xuu(i) 483 U = yuu(i) + dx*(Buu(i) + dx*(Cuu(i) + dx*Duu(i))) 484 endif 485 486 return 487end subroutine calc_U 488 489!------------------------------------------------------------------------------- 490! 491! Calculate embedding function U(rho) and first derivative dU(rho) 492! 493!------------------------------------------------------------------------------- 494! The "recursive" keyword is added below make this routine thread safe 495recursive subroutine calc_U_dU(rho,U,dU,ielast) 496 implicit none 497 498 !-- Transferred variables 499 real(c_double), intent(in) :: rho 500 real(c_double), intent(out) :: U,dU 501 integer(c_int), intent(inout) :: ielast 502 503 !-- Local variables 504 integer(c_int) i 505 real(c_double) dx 506 507 if (rho .le. rhomin) then 508 ! Argument less than the minimum stored value 509 U = 0.0_cd 510 dU = 0.0_cd 511 else 512 i = seval_i(nuu,rho,xuu,ielast) 513 dx = rho - xuu(i) 514 U = yuu(i) + dx*(Buu(i) + dx*(Cuu(i) + dx*Duu(i))) 515 dU = Buu(i) + dx*(2.0_cd*Cuu(i) + 3.0_cd*dx*Duu(i)) 516 endif 517 518 return 519end subroutine calc_U_dU 520 521!------------------------------------------------------------------------------- 522! 523! This function performs a binary search to find the index i 524! for evaluating the cubic spline function 525! 526! seval = y(i) + B(i)*(u-x(i)) + C(i)*(u-x(i))**2 + D(i)*(u-x(i))**3 527! 528! where x(i) .lt. u .lt. x(i+1), using horner's rule 529! 530! if u .lt. x(1) then i = 1 is used. 531! if u .ge. x(n) then i = n is used. 532! 533! input.. 534! 535! n = the number of data points 536! u = the abscissa at which the spline is to be evaluated 537! x = the array of data abscissas 538! i = current value of i 539! 540! if u is not in the same interval as the previous call, then a 541! binary search is performed to determine the proper interval. 542! 543!------------------------------------------------------------------------------- 544! The "recursive" keyword is added below make this routine thread safe 545recursive integer(c_int) function seval_i(n, u, x, i) 546 implicit none 547 548 !--Transferred variables 549 integer(c_int), intent(in) :: n 550 real(c_double), intent(in) :: u, x(n) 551 integer(c_int), intent(inout) :: i 552 553 !--Local variables 554 integer(c_int) j, k 555 556 if ( i .ge. n ) i = 1 557 if ( u .lt. x(i) ) go to 10 558 if ( u .le. x(i+1) ) go to 30 559 560 ! binary search 561 ! 562 10 i = 1 563 j = n+1 564 20 k = (i+j)/2 565 if ( u .lt. x(k) ) j = k 566 if ( u .ge. x(k) ) i = k 567 if ( j .gt. i+1 ) go to 20 568 569 ! got i, return 570 ! 571 30 seval_i = i 572 573 return 574 575end function seval_i 576 577!------------------------------------------------------------------------------- 578! 579! Compute energy and forces on atoms from the positions. 580! 581!------------------------------------------------------------------------------- 582! The "recursive" keyword is added below make this routine thread safe 583recursive subroutine Compute_Energy_Forces(model_compute_handle, & 584 model_compute_arguments_handle, ierr) bind(c) 585 implicit none 586 587 !-- Transferred variables 588 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle 589 type(kim_model_compute_arguments_handle_type), intent(in) :: & 590 model_compute_arguments_handle 591 integer(c_int), intent(out) :: ierr 592 593 !-- Local variables 594 real(c_double) :: Rij(DIM) 595 real(c_double) :: r,Rsqij,phi,dphi,g,dg,dU,U,dphieff 596 real(c_double) :: dphii,dUi,dphij,dUj 597 integer(c_int) :: i,j,jj,numnei,comp_force,comp_particleEnergy,comp_virial,comp_energy 598 integer(c_int) :: ierr2 599 real(c_double), allocatable :: rho(:),derU(:) 600 integer(c_int) :: irlast, ielast 601 602 !-- KIM variables 603 integer(c_int), pointer :: N 604 real(c_double), pointer :: energy 605 real(c_double), pointer :: coor(:,:) 606 real(c_double), pointer :: force(:,:) 607 real(c_double), pointer :: particleEnergy(:) 608 integer(c_int), pointer :: nei1part(:) 609 integer(c_int), pointer :: particleSpeciesCodes(:) 610 integer(c_int), pointer :: particleContributing(:) 611 real(c_double), pointer :: virial(:) 612 613 ! Avoid unused argument warning 614 if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue 615 616 ! Initialize interpolation interval indices 617 irlast = 1 618 ielast = 1 619 620 ! Unpack data from KIM object 621 ! 622 ierr = 0 623 call kim_get_argument_pointer( & 624 model_compute_arguments_handle, & 625 KIM_COMPUTE_ARGUMENT_NAME_NUMBER_OF_PARTICLES, N, ierr2) 626 ierr = ierr + ierr2 627 call kim_get_argument_pointer( & 628 model_compute_arguments_handle, & 629 KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_SPECIES_CODES, n, particleSpeciesCodes, & 630 ierr2) 631 ierr = ierr + ierr2 632 call kim_get_argument_pointer( & 633 model_compute_arguments_handle, & 634 KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_CONTRIBUTING, n, particleContributing, & 635 ierr2) 636 ierr = ierr + ierr2 637 call kim_get_argument_pointer( & 638 model_compute_arguments_handle, & 639 KIM_COMPUTE_ARGUMENT_NAME_COORDINATES, dim, n, coor, ierr2) 640 ierr = ierr + ierr2 641 call kim_get_argument_pointer( & 642 model_compute_arguments_handle, & 643 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, energy, ierr2) 644 ierr = ierr + ierr2 645 call kim_get_argument_pointer( & 646 model_compute_arguments_handle, & 647 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, dim, n, force, ierr2) 648 ierr = ierr + ierr2 649 call kim_get_argument_pointer( & 650 model_compute_arguments_handle, & 651 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_ENERGY, n, particleEnergy, ierr2) 652 ierr = ierr + ierr2 653 call kim_get_argument_pointer( & 654 model_compute_arguments_handle, & 655 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_VIRIAL, 6, virial, ierr2) 656 ierr = ierr + ierr2 657 if (ierr /= 0) then 658 call kim_log_entry(model_compute_arguments_handle, & 659 KIM_LOG_VERBOSITY_ERROR, & 660 "Failed to retrieve data from KIM API compute-arguments object") 661 return 662 endif 663 664 ! Check to see if we have been asked to compute the forces, energyperpart, 665 ! energy and virial 666 ! 667 if (associated(energy)) then 668 comp_energy = 1 669 else 670 comp_energy = 0 671 end if 672 if (associated(force)) then 673 comp_force = 1 674 else 675 comp_force = 0 676 end if 677 if (associated(particleEnergy)) then 678 comp_particleEnergy = 1 679 else 680 comp_particleEnergy = 0 681 end if 682 if (associated(virial)) then 683 comp_virial = 1 684 else 685 comp_virial = 0 686 end if 687 688 ! Check to be sure that the species are correct 689 ! 690 ierr = 1 ! assume an error 691 do i = 1,N 692 if (particleSpeciesCodes(i).ne.speccode) then 693 call kim_log_entry(model_compute_arguments_handle, & 694 KIM_LOG_VERBOSITY_ERROR, "Unexpected species code detected") 695 return 696 endif 697 enddo 698 ierr = 0 ! everything is ok 699 700 ! Initialize potential energies, forces, virial term, and electron density 701 ! 702 ! Note: that the variable `particleEnergy' does not need to be initialized 703 ! because it's initial value is set during the embedding energy 704 ! calculation. 705 if (comp_energy.eq.1) energy = 0.0_cd 706 if (comp_force.eq.1) force = 0.0_cd 707 if (comp_virial.eq.1) virial = 0.0_cd 708 dphieff = 0.0_cd 709 710 allocate( rho(N) ) ! pair functional electron density 711 rho = 0.0_cd 712 713 ! EAM embedded energy deriv 714 if (comp_force.eq.1.or.comp_virial.eq.1) allocate( derU(N) ) 715 716 ! Loop over particles in the neighbor list a first time, 717 ! to compute electron density (=coordination) 718 ! 719 do i = 1,N 720 if(particleContributing(i).eq.1) then 721 ! Get neighbor list of current atom 722 call kim_get_neighbor_list( & 723 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr) 724 if (ierr /= 0) then 725 ! some sort of problem, exit 726 call kim_log_entry(model_compute_arguments_handle, & 727 KIM_LOG_VERBOSITY_ERROR, "GetNeighborList failed") 728 ierr = 1 729 return 730 endif 731 732 ! Loop over the neighbors of atom i 733 ! 734 do jj = 1, numnei 735 736 j = nei1part(jj) ! get neighbor ID 737 if ( .not.( (particleContributing(j).eq.1) .and. & 738 j.lt.i) ) then ! effective half list 739 ! compute relative position vector 740 ! 741 Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j 742 743 ! compute contribution to electron density 744 ! 745 Rsqij = dot_product(Rij,Rij) ! compute square distance 746 if ( Rsqij .lt. model_cutsq ) then ! particles are interacting? 747 r = sqrt(Rsqij) ! compute distance 748 call calc_g(r,g,irlast) ! compute electron density 749 rho(i) = rho(i) + g ! accumulate electron density 750 if (particleContributing(j).eq.1) then 751 rho(j) = rho(j) + g ! this Model only supports a 752 ! single species, so we can 753 ! just add the same density 754 ! onto atom j 755 endif 756 endif 757 endif ! if ( i.lt.j ) 758 enddo ! loop on jj 759 endif ! Check on whether particle is contributing 760 enddo ! infinite do loop (terminated by exit statements above) 761 762 ! Now that we know the electron densities, calculate embedding part of energy 763 ! U and its derivative U' (derU) 764 ! 765 do i = 1,N 766 if(particleContributing(i).eq.1) then 767 if (comp_force.eq.1.or.comp_virial.eq.1) then 768 call calc_U_dU(rho(i),U,dU,ielast) ! compute embedding energy 769 ! and its derivative 770 derU(i) = dU ! store dU for later use 771 else 772 call calc_U(rho(i),U,ielast) ! compute just embedding energy 773 endif 774 775 ! accumulate the embedding energy contribution 776 ! 777 ! Assuming U(rho=0) = 0.0_cd 778 ! 779 if (comp_particleEnergy.eq.1) then ! accumulate embedding energy contribution 780 particleEnergy(i) = U 781 endif 782 if (comp_energy.eq.1) then 783 energy = energy + U 784 endif 785 786 endif ! Check on whether particle is contributing 787 enddo 788 789 ! Loop over particles in the neighbor list a second time, to compute 790 ! the forces and complete energy calculation 791 ! 792 793 do i = 1,N 794 if(particleContributing(i).eq.1) then 795 ! Get neighbor list of current atom 796 call kim_get_neighbor_list( & 797 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr) 798 if (ierr /= 0) then 799 ! some sort of problem, exit 800 call kim_log_entry(model_compute_arguments_handle, & 801 KIM_LOG_VERBOSITY_ERROR, "GetNeighborList failed") 802 ierr = 1 803 return 804 endif 805 806 ! Loop over the neighbors of atom i 807 ! 808 do jj = 1, numnei 809 810 j = nei1part(jj) ! get neighbor ID 811 if ( .not.( (particleContributing(j).eq.1) .and. & 812 j.lt.i) ) then ! effective half list 813 ! compute relative position vector 814 ! 815 Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j 816 817 ! compute energy and forces 818 ! 819 Rsqij = dot_product(Rij,Rij) ! compute square distance 820 if ( Rsqij .lt. model_cutsq ) then ! particles are interacting? 821 822 r = sqrt(Rsqij) ! compute distance 823 if (comp_force.eq.1.or.comp_virial.eq.1) then 824 call calc_phi_dphi(r,phi,dphi,irlast) ! compute pair potential 825 ! and it derivative 826 call calc_dg(r,dg,irlast) ! compute elect dens first deriv 827 828 if (particleContributing(j).eq.1) then 829 dphii = 0.5_cd*dphi 830 dphij = 0.5_cd*dphi 831 dUi = derU(i)*dg 832 dUj = derU(j)*dg 833 else 834 dphii = 0.5_cd*dphi 835 dphij = 0.0_cd 836 dUi = derU(i)*dg 837 dUj = 0.0_cd 838 endif 839 dphieff = dphii + dphij + dUi + dUj 840 else 841 call calc_phi(r,phi,irlast) ! compute just pair potential 842 endif 843 844 ! contribution to energy 845 ! 846 if (comp_particleEnergy.eq.1) then 847 particleEnergy(i) = particleEnergy(i) + 0.5_cd*phi ! accumulate energy 848 if (particleContributing(j).eq.1) then 849 particleEnergy(j) = particleEnergy(j) + 0.5_cd*phi ! accumulate energy 850 endif 851 endif 852 853 if (comp_energy.eq.1) then 854 if (particleContributing(j).eq.1) then 855 energy = energy + phi ! accumulate energy 856 else 857 energy = energy + 0.5_cd*phi 858 endif 859 endif 860 861 ! contribution to virial tensor 862 ! 863 if (comp_virial.eq.1) then 864 virial(1) = virial(1) + Rij(1)*Rij(1)*dphieff/r 865 virial(2) = virial(2) + Rij(2)*Rij(2)*dphieff/r 866 virial(3) = virial(3) + Rij(3)*Rij(3)*dphieff/r 867 virial(4) = virial(4) + Rij(2)*Rij(3)*dphieff/r 868 virial(5) = virial(5) + Rij(1)*Rij(3)*dphieff/r 869 virial(6) = virial(6) + Rij(1)*Rij(2)*dphieff/r 870 endif 871 872 ! contribution to forces 873 ! 874 if (comp_force.eq.1) then 875 force(:,i) = force(:,i) + dphieff*Rij/r ! accumulate force on atom i 876 force(:,j) = force(:,j) - dphieff*Rij/r ! accumulate force on atom j 877 endif 878 879 endif 880 endif ! if ( i.lt.j ) 881 enddo ! loop on jj 882 endif ! Check on whether particle is contributing 883 enddo ! infinite do loop (terminated by exit statements above) 884 885 ! Free temporary storage 886 ! 887 deallocate( rho ) 888 if (comp_force.eq.1.or.comp_virial.eq.1) deallocate( derU ) 889 890 ! Everything is great 891 ierr = 0 892 893 return 894end subroutine Compute_Energy_Forces 895 896!------------------------------------------------------------------------------- 897! 898! Model destroy routine 899! 900!------------------------------------------------------------------------------- 901! The "recursive" keyword is added below make this routine thread safe 902recursive subroutine destroy(model_destroy_handle, ierr) bind(c) 903 use, intrinsic :: iso_c_binding 904 implicit none 905 906 !-- Transferred variables 907 type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle 908 integer(c_int), intent(out) :: ierr 909 type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf 910 911 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf) 912 call c_f_pointer(pbuf, buf) 913 914 call kim_log_entry(model_destroy_handle, KIM_LOG_VERBOSITY_INFORMATION, & 915 "deallocating model buffer") 916 deallocate(buf) 917 ierr = 0 ! everything is good 918 919 return 920end subroutine destroy 921 922!------------------------------------------------------------------------------- 923! 924! Model compute arguments create routine (REQUIRED) 925! 926!------------------------------------------------------------------------------- 927! The "recursive" keyword is added below make this routine thread safe 928recursive subroutine compute_arguments_create(model_compute_handle, & 929 model_compute_arguments_create_handle, ierr) bind(c) 930 use, intrinsic :: iso_c_binding 931 use kim_model_compute_arguments_create_module 932 implicit none 933 934 !-- Transferred variables 935 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle 936 type(kim_model_compute_arguments_create_handle_type), intent(inout) :: & 937 model_compute_arguments_create_handle 938 integer(c_int), intent(out) :: ierr 939 940 integer(c_int) :: ierr2 941 942 ! Avoid unused argument warning 943 if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue 944 945 ierr = 0 946 ierr2 = 0 947 948 ! register arguments 949 call kim_set_argument_support_status( & 950 model_compute_arguments_create_handle, & 951 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, & 952 KIM_SUPPORT_STATUS_OPTIONAL, ierr2) 953 ierr = ierr + ierr2 954 call kim_set_argument_support_status( & 955 model_compute_arguments_create_handle, & 956 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, & 957 KIM_SUPPORT_STATUS_OPTIONAL, ierr2) 958 ierr = ierr + ierr2 959 call kim_set_argument_support_status( & 960 model_compute_arguments_create_handle, & 961 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_ENERGY, & 962 KIM_SUPPORT_STATUS_OPTIONAL, ierr2) 963 ierr = ierr + ierr2 964 call kim_set_argument_support_status( & 965 model_compute_arguments_create_handle, & 966 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_VIRIAL, & 967 KIM_SUPPORT_STATUS_OPTIONAL, ierr2) 968 ierr = ierr + ierr2 969 970 ! register call backs (ProcessDEDrTerm, ProcessD2EDr2Term 971 ! NONE 972 973 if (ierr /= 0) then 974 ierr = 1 975 call kim_log_entry(model_compute_arguments_create_handle, & 976 KIM_LOG_VERBOSITY_ERROR, & 977 "Unable to successfully create compute-arguments object") 978 endif 979 980 return 981end subroutine compute_arguments_create 982 983!------------------------------------------------------------------------------- 984! 985! Model compute arguments destroy routine (REQUIRED) 986! 987!------------------------------------------------------------------------------- 988! The "recursive" keyword is added below make this routine thread safe 989recursive subroutine compute_arguments_destroy(model_compute_handle, & 990 model_compute_arguments_destroy_handle, ierr) bind(c) 991 use, intrinsic :: iso_c_binding 992 implicit none 993 994 !-- Transferred variables 995 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle 996 type(kim_model_compute_arguments_destroy_handle_type), intent(inout) :: & 997 model_compute_arguments_destroy_handle 998 integer(c_int), intent(out) :: ierr 999 1000 ! avoid unused dummy argument warnings 1001 if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue 1002 if (model_compute_arguments_destroy_handle .eq. & 1003 KIM_MODEL_COMPUTE_ARGUMENTS_DESTROY_NULL_HANDLE) continue 1004 1005 ierr = 0 1006 return 1007end subroutine compute_arguments_destroy 1008 1009end module EAM_ErcolessiAdams_1994_Al 1010 1011!------------------------------------------------------------------------------- 1012! 1013! Model initialization routine (REQUIRED) 1014! 1015!------------------------------------------------------------------------------- 1016! The "recursive" keyword is added below make this routine thread safe 1017recursive subroutine create(model_create_handle, requested_length_unit, & 1018 requested_energy_unit, requested_charge_unit, requested_temperature_unit, & 1019 requested_time_unit, ierr) bind(c) 1020 use, intrinsic :: iso_c_binding 1021 use EAM_ErcolessiAdams_1994_Al 1022 use kim_model_headers_module 1023 implicit none 1024 1025 !-- Transferred variables 1026 type(kim_model_create_handle_type), intent(inout) :: model_create_handle 1027 type(kim_length_unit_type), intent(in) :: requested_length_unit 1028 type(kim_energy_unit_type), intent(in) :: requested_energy_unit 1029 type(kim_charge_unit_type), intent(in) :: requested_charge_unit 1030 type(kim_temperature_unit_type), intent(in) :: requested_temperature_unit 1031 type(kim_time_unit_type), intent(in) :: requested_time_unit 1032 integer(c_int), intent(out) :: ierr 1033 1034 !-- KIM variables 1035 integer(c_int) :: ierr2 1036 type(buffer_type), pointer :: buf 1037 1038 ierr = 0 1039 ierr2 = 0 1040 1041 call kim_set_units(model_create_handle, & 1042 KIM_LENGTH_UNIT_A, & 1043 KIM_ENERGY_UNIT_EV, & 1044 KIM_CHARGE_UNIT_UNUSED, & 1045 KIM_TEMPERATURE_UNIT_UNUSED, & 1046 KIM_TIME_UNIT_UNUSED, & 1047 ierr2) 1048 ierr = ierr + ierr2 1049 1050 ! register species 1051 call kim_set_species_code(model_create_handle, & 1052 KIM_SPECIES_NAME_AL, speccode, ierr2) 1053 ierr = ierr + ierr2 1054 1055 ! register numbering 1056 call kim_set_model_numbering(model_create_handle, & 1057 KIM_NUMBERING_ONE_BASED, ierr2); 1058 ierr = ierr + ierr2 1059 1060 ! register function pointers 1061 call kim_set_routine_pointer(model_create_handle, & 1062 KIM_MODEL_ROUTINE_NAME_COMPUTE, KIM_LANGUAGE_NAME_FORTRAN, & 1063 1, c_funloc(Compute_Energy_Forces), ierr2) 1064 ierr = ierr + ierr2 1065 call kim_set_routine_pointer( & 1066 model_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_CREATE, & 1067 KIM_LANGUAGE_NAME_FORTRAN, 1, c_funloc(compute_arguments_create), ierr2) 1068 ierr = ierr + ierr2 1069 call kim_set_routine_pointer( & 1070 model_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_DESTROY, & 1071 KIM_LANGUAGE_NAME_FORTRAN, 1, c_funloc(compute_arguments_destroy), ierr2) 1072 ierr = ierr + ierr2 1073 call kim_set_routine_pointer(model_create_handle, & 1074 KIM_MODEL_ROUTINE_NAME_DESTROY, KIM_LANGUAGE_NAME_FORTRAN, & 1075 1, c_funloc(destroy), ierr2) 1076 ierr = ierr + ierr2 1077 1078 ! allocate buffer 1079 allocate(buf) 1080 1081 ! store model buffer in KIM object 1082 call kim_set_model_buffer_pointer(model_create_handle, & 1083 c_loc(buf)) 1084 1085 ! set buffer values 1086 buf%influence_distance = model_cutoff 1087 buf%cutoff(1) = model_cutoff 1088 buf%model_will_not_request_neighbors_of_noncontributing_particles(1) = 1 1089 1090 ! register influence distance 1091 call kim_set_influence_distance_pointer( & 1092 model_create_handle, buf%influence_distance) 1093 1094 ! register cutoff 1095 call kim_set_neighbor_list_pointers(model_create_handle, & 1096 1, buf%cutoff, & 1097 buf%model_will_not_request_neighbors_of_noncontributing_particles) 1098 1099 if (ierr /= 0) then 1100 ierr = 1 1101 call kim_log_entry(model_create_handle, KIM_LOG_VERBOSITY_ERROR, & 1102 "Unable to successfully initialize model") 1103 endif 1104 1105 return 1106end subroutine create 1107