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! Ellad B. Tadmor 27! Ryan S. Elliott 28! Valeriu Smirichinski 29! Stephen M. Whalen 30! 31 32!**************************************************************************** 33!** 34!** MODULE Pair_Lennard_Jones_Shifted 35!** 36!** Lennard-Jones pair potential KIM Model Driver 37!** shifted to have zero energy at the cutoff radius 38!** 39!** Language: Fortran 2003 40!** 41!**************************************************************************** 42module Pair_Lennard_Jones_Shifted 43 44use, intrinsic :: iso_c_binding 45use kim_model_driver_headers_module 46implicit none 47 48save 49private 50public BUFFER_TYPE, & 51 compute, & 52 compute_arguments_create, & 53 compute_arguments_destroy, & 54 refresh, & 55 destroy, & 56 calc_phi, & 57 calc_phi_dphi, & 58 calc_phi_dphi_d2phi, & 59 speccode 60 61! Below are the definitions and values of all Model parameters 62integer(c_int), parameter :: cd = c_double ! for literal constants 63integer(c_int), parameter :: DIM=3 ! dimensionality of space 64integer(c_int), parameter :: speccode = 1 ! internal species code 65 66!------------------------------------------------------------------------------- 67! 68! Definition of Buffer type 69! 70!------------------------------------------------------------------------------- 71type, bind(c) :: BUFFER_TYPE 72 real(c_double) :: influence_distance 73 real(c_double) :: cutoff(1) 74 real(c_double) :: cutsq(1) 75 real(c_double) :: eps(1) 76 real(c_double) :: sigma(1) 77 real(c_double) :: shift(1) 78 integer(c_int) :: & 79 model_will_not_request_neighbors_of_noncontributing_particles(1) 80endtype BUFFER_TYPE 81 82 83contains 84 85!------------------------------------------------------------------------------- 86! 87! Calculate pair potential phi(r) 88! 89!------------------------------------------------------------------------------- 90! The "recursive" keyword is added below make this routine thread safe 91recursive subroutine calc_phi(model_eps, & 92 model_sigma, & 93 model_shift, & 94 model_cutoff,r,phi) 95 implicit none 96 97 !-- Transferred variables 98 real(c_double), intent(in) :: model_eps 99 real(c_double), intent(in) :: model_sigma 100 real(c_double), intent(in) :: model_shift 101 real(c_double), intent(in) :: model_cutoff 102 real(c_double), intent(in) :: r 103 real(c_double), intent(out) :: phi 104 105 !-- Local variables 106 real(c_double) rsq,sor,sor6,sor12 107 108 if (r .gt. model_cutoff) then 109 ! Argument exceeds cutoff radius 110 phi = 0.0_cd 111 return 112 endif 113 114 rsq = r*r ! r^2 115 sor = model_sigma/r ! (sig/r) 116 sor6 = sor*sor*sor ! 117 sor6 = sor6*sor6 ! (sig/r)^6 118 sor12 = sor6*sor6 ! (sig/r)^12 119 phi = 4.0_cd*model_eps*(sor12-sor6) + model_shift 120 121 return 122end subroutine calc_phi 123 124!------------------------------------------------------------------------------- 125! 126! Calculate pair potential phi(r) and its derivative dphi(r) 127! 128!------------------------------------------------------------------------------- 129! The "recursive" keyword is added below make this routine thread safe 130recursive subroutine calc_phi_dphi(model_eps, & 131 model_sigma, & 132 model_shift, & 133 model_cutoff,r,phi,dphi) 134 implicit none 135 136 !-- Transferred variables 137 real(c_double), intent(in) :: model_eps 138 real(c_double), intent(in) :: model_sigma 139 real(c_double), intent(in) :: model_shift 140 real(c_double), intent(in) :: model_cutoff 141 real(c_double), intent(in) :: r 142 real(c_double), intent(out) :: phi,dphi 143 144 !-- Local variables 145 real(c_double) rsq,sor,sor6,sor12 146 147 if (r .gt. model_cutoff) then 148 ! Argument exceeds cutoff radius 149 phi = 0.0_cd 150 dphi = 0.0_cd 151 return 152 endif 153 154 rsq = r*r ! r^2 155 sor = model_sigma/r ! (sig/r) 156 sor6 = sor*sor*sor ! 157 sor6 = sor6*sor6 ! (sig/r)^6 158 sor12= sor6*sor6 ! (sig/r)^12 159 phi = 4.0_cd*model_eps*(sor12-sor6) + model_shift 160 dphi = 24.0_cd*model_eps*(-2.0_cd*sor12+sor6)/r 161 162 return 163end subroutine calc_phi_dphi 164 165!------------------------------------------------------------------------------- 166! 167! Calculate pair potential phi(r) and its derivatives dphi(r) and d2phi(r) 168! 169!------------------------------------------------------------------------------- 170! The "recursive" keyword is added below make this routine thread safe 171recursive subroutine calc_phi_dphi_d2phi(model_eps, & 172 model_sigma, & 173 model_shift, & 174 model_cutoff,r,phi,dphi,d2phi) 175 implicit none 176 177 !-- Transferred variables 178 real(c_double), intent(in) :: model_eps 179 real(c_double), intent(in) :: model_sigma 180 real(c_double), intent(in) :: model_shift 181 real(c_double), intent(in) :: model_cutoff 182 real(c_double), intent(in) :: r 183 real(c_double), intent(out) :: phi,dphi,d2phi 184 185 !-- Local variables 186 real(c_double) rsq,sor,sor6,sor12 187 188 if (r .gt. model_cutoff) then 189 ! Argument exceeds cutoff radius 190 phi = 0.0_cd 191 dphi = 0.0_cd 192 d2phi = 0.0_cd 193 return 194 endif 195 196 rsq = r*r ! r^2 197 sor = model_sigma/r ! (sig/r) 198 sor6 = sor*sor*sor ! 199 sor6 = sor6*sor6 ! (sig/r)^6 200 sor12 = sor6*sor6 ! (sig/r)^12 201 phi = 4.0_cd*model_eps*(sor12-sor6) + model_shift 202 dphi = 24.0_cd*model_eps*(-2.0_cd*sor12+sor6)/r 203 d2phi = 24.0_cd*model_eps*(26.0_cd*sor12-7.0_cd*sor6)/rsq 204 205 return 206end subroutine calc_phi_dphi_d2phi 207 208!------------------------------------------------------------------------------- 209! 210! Compute 211! 212!------------------------------------------------------------------------------- 213! The "recursive" keyword is added below make this routine thread safe 214recursive subroutine compute(model_compute_handle, & 215 model_compute_arguments_handle, ierr) bind(c) 216 implicit none 217 218 !-- Transferred variables 219 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle 220 type(kim_model_compute_arguments_handle_type), intent(in) :: & 221 model_compute_arguments_handle 222 integer(c_int), intent(out) :: ierr 223 224 !-- Local variables 225 real(c_double) :: r,Rsqij,phi,dphi,d2phi,dEidr,d2Eidr 226 real(c_double) :: v 227 real(c_double) :: vir(6) 228 integer(c_int) :: i,j,jj,numnei 229 integer(c_int) :: ierr2 230 integer(c_int) :: comp_force,comp_energy,comp_particleEnergy,comp_process_dEdr, & 231 comp_process_d2Edr2,comp_virial,comp_particleVirial 232 type(BUFFER_TYPE), pointer :: buf; type(c_ptr) :: pbuf 233 234 real(c_double), target :: Rij(DIM) 235 ! Quantities for computing d2Edr2 236 real(c_double), target :: Rij_pairs(DIM,2) 237 real(c_double), target :: r_pairs(2) 238 integer(c_int), target :: i_pairs(2), j_pairs(2) 239 240 !-- KIM variables 241 integer(c_int), pointer :: N 242 real(c_double), pointer :: energy 243 real(c_double), pointer :: coor(:,:) 244 real(c_double), pointer :: force(:,:) 245 real(c_double), pointer :: particleEnergy(:) 246 real(c_double), pointer :: virial(:) 247 real(c_double), pointer :: particleVirial(:, :) 248 integer(c_int), pointer :: nei1part(:) 249 integer(c_int), pointer :: particleSpeciesCodes(:) 250 integer(c_int), pointer :: particleContributing(:) 251 252 ! get model buffer from KIM object 253 call kim_get_model_buffer_pointer(model_compute_handle, pbuf) 254 call c_f_pointer(pbuf, buf) 255 256 ! Unpack data from KIM object 257 ! 258 ierr = 0_c_int 259 call kim_get_argument_pointer( & 260 model_compute_arguments_handle, & 261 KIM_COMPUTE_ARGUMENT_NAME_NUMBER_OF_PARTICLES, N, ierr2) 262 ierr = ierr + ierr2 263 264 call kim_get_argument_pointer( & 265 model_compute_arguments_handle, & 266 KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_SPECIES_CODES, & 267 N, particlespeciesCodes, ierr2) 268 ierr = ierr + ierr2 269 call kim_get_argument_pointer( & 270 model_compute_arguments_handle, & 271 KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_CONTRIBUTING, N, particlecontributing, & 272 ierr2) 273 ierr = ierr + ierr2 274 call kim_get_argument_pointer( & 275 model_compute_arguments_handle, & 276 KIM_COMPUTE_ARGUMENT_NAME_COORDINATES, DIM, N, coor, ierr2) 277 ierr = ierr + ierr2 278 call kim_get_argument_pointer( & 279 model_compute_arguments_handle, & 280 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, energy, ierr2) 281 ierr = ierr + ierr2 282 call kim_get_argument_pointer( & 283 model_compute_arguments_handle, & 284 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, DIM, N, force, ierr2) 285 ierr = ierr + ierr2 286 call kim_get_argument_pointer( & 287 model_compute_arguments_handle, & 288 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_ENERGY, N, particleEnergy, ierr2) 289 ierr = ierr + ierr2 290 291 call kim_get_argument_pointer( & 292 model_compute_arguments_handle, & 293 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_VIRIAL, & 294 6_c_int, & 295 virial, & 296 ierr2) 297 ierr = ierr + ierr2 298 299 call kim_get_argument_pointer( & 300 model_compute_arguments_handle, & 301 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_VIRIAL, & 302 6_c_int, & 303 N, & 304 particleVirial, & 305 ierr2) 306 ierr = ierr + ierr2 307 308 if (ierr /= 0_c_int) then 309 call kim_log_entry(model_compute_arguments_handle, & 310 KIM_LOG_VERBOSITY_ERROR, "get_argument_pointer") 311 return 312 endif 313 314 ! Check to see if we have been asked to compute the energy, forces, energy per particle, 315 ! dEdr, and d2Edr2 316 ! 317 if (associated(energy)) then 318 comp_energy = 1_c_int 319 else 320 comp_energy = 0_c_int 321 end if 322 if (associated(force)) then 323 comp_force = 1_c_int 324 else 325 comp_force = 0_c_int 326 end if 327 if (associated(particleEnergy)) then 328 comp_particleEnergy = 1_c_int 329 else 330 comp_particleEnergy = 0_c_int 331 end if 332 if (associated(virial)) then 333 comp_virial = 1_c_int 334 else 335 comp_virial = 0_c_int 336 end if 337 if (associated(particleVirial)) then 338 comp_particleVirial = 1_c_int 339 else 340 comp_particleVirial = 0_c_int 341 end if 342 343 ierr = 0_c_int 344 call kim_is_callback_present( & 345 model_compute_arguments_handle, & 346 KIM_COMPUTE_CALLBACK_NAME_PROCESS_DEDR_TERM, comp_process_dedr, ierr2) 347 ierr = ierr + ierr2 348 call kim_is_callback_present( & 349 model_compute_arguments_handle, & 350 KIM_COMPUTE_CALLBACK_NAME_PROCESS_D2edr2_term, comp_process_d2edr2, ierr2) 351 ierr = ierr + ierr2 352 if (ierr /= 0_c_int) then 353 call kim_log_entry(model_compute_arguments_handle, & 354 KIM_LOG_VERBOSITY_ERROR, "get_compute") 355 return 356 endif 357 358 ! Check to be sure that the species are correct 359 ! 360 ierr = 1_c_int ! assume an error 361 do i = 1_c_int,N 362 if (particleSpeciesCodes(i).ne.speccode) then 363 call kim_log_entry(model_compute_arguments_handle, & 364 KIM_LOG_VERBOSITY_ERROR, "Unexpected species code detected") 365 return 366 endif 367 enddo 368 ierr = 0_c_int ! everything is ok 369 370 ! Initialize potential energies, forces 371 ! 372 if (comp_particleEnergy.eq.1_c_int) particleEnergy = 0.0_cd 373 if (comp_energy.eq.1_c_int) energy = 0.0_cd 374 if (comp_force.eq.1_c_int) force = 0.0_cd 375 if (comp_virial.eq.1_c_int) virial = 0.0_cd 376 if (comp_particleVirial.eq.1_c_int) particleVirial = 0.0_cd 377 378 ! 379 ! Compute energy and forces 380 ! 381 382 ! Loop over particles and compute energy and forces 383 ! 384 do i = 1_c_int,N 385 if (particleContributing(i) == 1_c_int) then 386 ! Set up neighbor list for next particle 387 call kim_get_neighbor_list( & 388 model_compute_arguments_handle, 1_c_int, i, numnei, nei1part, ierr) 389 if (ierr /= 0_c_int) then 390 ! some sort of problem, exit 391 call kim_log_entry(model_compute_arguments_handle, & 392 KIM_LOG_VERBOSITY_ERROR, "GetNeighborList failed") 393 ierr = 1_c_int 394 return 395 endif 396 397 ! Loop over the neighbors of atom i 398 ! 399 do jj = 1_c_int, numnei 400 401 j = nei1part(jj) ! get neighbor ID 402 403 if ( .not.( (particleContributing(j) == 1_c_int) .and. & 404 j.lt.i) ) then ! effective half list 405 ! compute relative position vector 406 ! 407 Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j 408 409 ! compute energy and forces 410 ! 411 Rsqij = dot_product(Rij,Rij) ! compute square distance 412 if ( Rsqij .lt. buf%cutsq(1) ) then ! particles are interacting? 413 r = sqrt(Rsqij) ! compute distance 414 if (comp_process_d2Edr2.eq.1_c_int) then 415 call calc_phi_dphi_d2phi(buf%eps(1), & 416 buf%sigma(1), & 417 buf%shift(1), & 418 buf%cutoff(1), & 419 r,phi,dphi,d2phi) ! compute pair potential 420 ! and its derivatives 421 if (particleContributing(j).eq.1_c_int) then 422 dEidr = dphi 423 d2Eidr = d2phi 424 else 425 dEidr = 0.5_cd*dphi 426 d2Eidr = 0.5_cd*d2phi 427 endif 428 elseif ((comp_force == 1_c_int) .or. & 429 (comp_process_dEdr == 1_c_int) .or. & 430 (comp_virial == 1_c_int) .or. & 431 (comp_particleVirial== 1_c_int)) then 432 call calc_phi_dphi(buf%eps(1), & 433 buf%sigma(1), & 434 buf%shift(1), & 435 buf%cutoff(1), & 436 r,phi,dphi) ! compute pair potential 437 ! and it derivative 438 439 if (particleContributing(j).eq.1_c_int) then 440 dEidr = dphi 441 else 442 dEidr = 0.5_cd*dphi 443 endif 444 else 445 call calc_phi(buf%eps(1), & 446 buf%sigma(1), & 447 buf%shift(1), & 448 buf%cutoff(1), & 449 r,phi) ! compute just pair potential 450 endif 451 452 ! contribution to energy 453 ! 454 if (comp_particleEnergy.eq.1_c_int) then 455 particleEnergy(i) = particleEnergy(i) + 0.5_cd*phi ! accumulate energy 456 if (particleContributing(j).eq.1_c_int) then 457 particleEnergy(j) = particleEnergy(j) + 0.5_cd*phi ! accumulate energy 458 endif 459 endif 460 461 if (comp_energy.eq.1_c_int) then 462 if (particleContributing(j).eq.1_c_int) then 463 energy = energy + phi ! add v to total energy 464 else 465 energy = energy + 0.5_cd*phi ! Throw out the ghost atom half of the energy 466 endif 467 endif 468 469 ! contribution to process_dEdr 470 ! 471 if (comp_process_dEdr.eq.1_c_int) then 472 call kim_process_dedr_term( & 473 model_compute_arguments_handle, dEidr, r, Rij, i, j, ierr) 474 endif 475 476 ! contribution to virial and particle virial 477 ! 478 if ((comp_virial == 1_c_int) .or. & 479 (comp_particleVirial == 1_c_int)) then 480 ! Virial has 6 components and is stored as a 6-element 481 ! vector in the following order: xx, yy, zz, yz, xz, xy. 482 v = dEidr / r 483 484 vir(1) = v * Rij(1) * Rij(1) 485 vir(2) = v * Rij(2) * Rij(2) 486 vir(3) = v * Rij(3) * Rij(3) 487 vir(4) = v * Rij(2) * Rij(3) 488 vir(5) = v * Rij(1) * Rij(3) 489 vir(6) = v * Rij(1) * Rij(2) 490 491 if (comp_virial == 1_c_int) then 492 virial(1) = virial(1) + vir(1) 493 virial(2) = virial(2) + vir(2) 494 virial(3) = virial(3) + vir(3) 495 virial(4) = virial(4) + vir(4) 496 virial(5) = virial(5) + vir(5) 497 virial(6) = virial(6) + vir(6) 498 endif 499 500 if (comp_particleVirial == 1_c_int) then 501 vir = vir * 0.5 502 503 particleVirial(1, i) = particleVirial(1, i) + vir(1) 504 particleVirial(2, i) = particleVirial(2, i) + vir(2) 505 particleVirial(3, i) = particleVirial(3, i) + vir(3) 506 particleVirial(4, i) = particleVirial(4, i) + vir(4) 507 particleVirial(5, i) = particleVirial(5, i) + vir(5) 508 particleVirial(6, i) = particleVirial(6, i) + vir(6) 509 510 particleVirial(1, j) = particleVirial(1, j) + vir(1) 511 particleVirial(2, j) = particleVirial(2, j) + vir(2) 512 particleVirial(3, j) = particleVirial(3, j) + vir(3) 513 particleVirial(4, j) = particleVirial(4, j) + vir(4) 514 particleVirial(5, j) = particleVirial(5, j) + vir(5) 515 particleVirial(6, j) = particleVirial(6, j) + vir(6) 516 endif 517 endif 518 519 ! contribution to process_d2Edr2 520 if (comp_process_d2Edr2.eq.1_c_int) then 521 r_pairs(1) = r 522 r_pairs(2) = r 523 Rij_pairs(:,1) = Rij 524 Rij_pairs(:,2) = Rij 525 i_pairs(1) = i 526 i_pairs(2) = i 527 j_pairs(1) = j 528 j_pairs(2) = j 529 530 call kim_process_d2edr2_term( & 531 model_compute_arguments_handle, d2Eidr, & 532 r_pairs, Rij_pairs, i_pairs, j_pairs, ierr) 533 534 endif 535 536 ! contribution to forces 537 ! 538 if (comp_force.eq.1_c_int) then 539 force(:,i) = force(:,i) + dEidr*Rij/r ! accumulate force on atom i 540 force(:,j) = force(:,j) - dEidr*Rij/r ! accumulate force on atom j 541 endif 542 543 endif ! Check on whether particle jj is interacting with particle i 544 endif ! if ( i.lt.j ) 545 enddo ! loop on jj 546 endif ! Check on whether particle i is contributing 547 enddo ! infinite do loop (terminated by exit statements above) 548 549 ! No errors 550 ! 551 ierr = 0_c_int 552 return 553 554end subroutine compute 555 556!------------------------------------------------------------------------------- 557! 558! Model driver refresh routine 559! 560!------------------------------------------------------------------------------- 561! The "recursive" keyword is added below make this routine thread safe 562recursive subroutine refresh(model_refresh_handle, ierr) bind(c) 563 implicit none 564 565 !-- Transferred variables 566 type(kim_model_refresh_handle_type), intent(inout) :: model_refresh_handle 567 integer(c_int), intent(out) :: ierr 568 569 !-- Local variables 570 real(c_double) energy_at_cutoff 571 type(BUFFER_TYPE), pointer :: buf; type(c_ptr) :: pbuf 572 573 ! Get model buffer from KIM object 574 call kim_get_model_buffer_pointer(model_refresh_handle, pbuf) 575 call c_f_pointer(pbuf, buf) 576 577 call kim_set_influence_distance_pointer(model_refresh_handle, & 578 buf%influence_distance) 579 call kim_set_neighbor_list_pointers(model_refresh_handle, & 580 1_c_int, buf%cutoff, & 581 buf%model_will_not_request_neighbors_of_noncontributing_particles) 582 583 ! Set new values in KIM object and buffer 584 ! 585 buf%influence_distance = buf%cutoff(1) 586 buf%cutsq(1) = buf%cutoff(1)*buf%cutoff(1) 587 588 ! Calculate pair potential at r=cutoff with shift=0.0 589 call calc_phi(buf%eps(1), & 590 buf%sigma(1), & 591 0.0_cd, & 592 buf%cutoff(1), & 593 buf%cutoff(1),energy_at_cutoff) 594 595 ! Set the corresponding necessary shift into the 'shift' buffer variable 596 buf%shift(1) = -energy_at_cutoff 597 598 ierr = 0_c_int 599 return 600 601end subroutine refresh 602 603!------------------------------------------------------------------------------- 604! 605! Model driver destroy routine 606! 607!------------------------------------------------------------------------------- 608! The "recursive" keyword is added below make this routine thread safe 609recursive subroutine destroy(model_destroy_handle, ierr) bind(c) 610 implicit none 611 612 !-- Transferred variables 613 type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle 614 integer(c_int), intent(out) :: ierr 615 616 !-- Local variables 617 type(BUFFER_TYPE), pointer :: buf; type(c_ptr) :: pbuf 618 619 ! get model buffer from KIM object 620 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf) 621 call c_f_pointer(pbuf, buf) 622 623 deallocate( buf ) 624 625 ierr = 0_c_int 626 return 627 628end subroutine destroy 629 630!------------------------------------------------------------------------------- 631! 632! Model driver compute arguments create routine 633! 634!------------------------------------------------------------------------------- 635! The "recursive" keyword is added below make this routine thread safe 636recursive subroutine compute_arguments_create(model_compute_handle, & 637 model_compute_arguments_create_handle, ierr) bind(c) 638use kim_model_compute_arguments_create_module 639 implicit none 640 641 !-- Transferred variables 642 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle 643 type(kim_model_compute_arguments_create_handle_type), intent(inout) :: & 644 model_compute_arguments_create_handle 645 integer(c_int), intent(out) :: ierr 646 647 integer(c_int) ierr2 648 649 ! avoid unused dummy argument warnings 650 if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue 651 652 ierr = 0_c_int 653 ierr2 = 0_c_int 654 655 ! register arguments 656 call kim_set_argument_support_status( & 657 model_compute_arguments_create_handle, & 658 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, & 659 KIM_SUPPORT_STATUS_OPTIONAL, ierr) 660 call kim_set_argument_support_status( & 661 model_compute_arguments_create_handle, & 662 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, & 663 KIM_SUPPORT_STATUS_OPTIONAL, ierr2) 664 ierr = ierr + ierr2 665 call kim_set_argument_support_status( & 666 model_compute_arguments_create_handle, & 667 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_ENERGY, & 668 KIM_SUPPORT_STATUS_OPTIONAL, ierr2) 669 ierr = ierr + ierr2 670 call kim_set_argument_support_status( & 671 model_compute_arguments_create_handle, & 672 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_VIRIAL, & 673 KIM_SUPPORT_STATUS_OPTIONAL, & 674 ierr2) 675 ierr = ierr + ierr2 676 call kim_set_argument_support_status( & 677 model_compute_arguments_create_handle, & 678 KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_VIRIAL, & 679 KIM_SUPPORT_STATUS_OPTIONAL, & 680 ierr2) 681 ierr = ierr + ierr2 682 if (ierr /= 0_c_int) then 683 call kim_log_entry(model_compute_arguments_create_handle, & 684 KIM_LOG_VERBOSITY_ERROR, "Unable to register arguments support status") 685 goto 42 686 end if 687 688 ! register callbacks 689 call kim_set_callback_support_status( & 690 model_compute_arguments_create_handle, & 691 KIM_COMPUTE_CALLBACK_NAME_PROCESS_DEDR_TERM, & 692 KIM_SUPPORT_STATUS_OPTIONAL, ierr) 693 call kim_set_callback_support_status( & 694 model_compute_arguments_create_handle, & 695 KIM_COMPUTE_CALLBACK_NAME_PROCESS_D2edr2_term, & 696 KIM_SUPPORT_STATUS_OPTIONAL, ierr2) 697 ierr = ierr + ierr2 698 if (ierr /= 0_c_int) then 699 call kim_log_entry(model_compute_arguments_create_handle, & 700 KIM_LOG_VERBOSITY_ERROR, "Unable to register callback support status") 701 goto 42 702 end if 703 704 ierr = 0_c_int 705 42 continue 706 return 707 708end subroutine compute_arguments_create 709 710!------------------------------------------------------------------------------- 711! 712! Model driver compute arguments destroy routine 713! 714!------------------------------------------------------------------------------- 715! The "recursive" keyword is added below make this routine thread safe 716recursive subroutine compute_arguments_destroy(model_compute_handle, & 717 model_compute_arguments_destroy_handle, ierr) bind(c) 718 implicit none 719 720 !-- Transferred variables 721 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle 722 type(kim_model_compute_arguments_destroy_handle_type), intent(in) :: & 723 model_compute_arguments_destroy_handle 724 integer(c_int), intent(out) :: ierr 725 726 ! avoid unsed dummy argument warnings 727 if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue 728 if (model_compute_arguments_destroy_handle .eq. & 729 KIM_MODEL_COMPUTE_ARGUMENTS_DESTROY_NULL_HANDLE) continue 730 731 ierr = 0_c_int 732 return 733 734end subroutine compute_arguments_destroy 735 736end module Pair_Lennard_Jones_Shifted 737 738!------------------------------------------------------------------------------- 739! 740! Model driver initialization routine (REQUIRED) 741! 742!------------------------------------------------------------------------------- 743! The "recursive" keyword is added below make this routine thread safe 744recursive subroutine create(model_driver_create_handle, requested_length_unit, & 745 requested_energy_unit, requested_charge_unit, requested_temperature_unit, & 746 requested_time_unit, ierr) bind(c) 747 use, intrinsic :: iso_c_binding 748 use Pair_Lennard_Jones_Shifted 749 use kim_model_driver_headers_module 750 implicit none 751 752 !-- Transferred variables 753 type(kim_model_driver_create_handle_type), intent(inout) :: model_driver_create_handle 754 type(kim_length_unit_type), intent(in), value :: requested_length_unit 755 type(kim_energy_unit_type), intent(in), value :: requested_energy_unit 756 type(kim_charge_unit_type), intent(in), value :: requested_charge_unit 757 type(kim_temperature_unit_type), intent(in), value :: requested_temperature_unit 758 type(kim_time_unit_type), intent(in), value :: requested_time_unit 759 integer(c_int), intent(out) :: ierr 760 761 !-- Local variables 762 integer(c_int), parameter :: cd = c_double ! used for literal constants 763 integer(c_int) :: number_of_parameter_files 764 character(len=1024, kind=c_char) :: parameter_file_name 765 integer(c_int) :: ierr2 766 character(len=100, kind=c_char) :: in_species 767 real(c_double) :: in_cutoff 768 real(c_double) :: in_eps 769 real(c_double) :: in_sigma 770 real(c_double) :: energy_at_cutoff 771 real(c_double) :: convert_length, convert_energy 772 type(kim_species_name_type) species_name 773 type(BUFFER_TYPE), pointer :: buf 774 775 ! Register numbering 776 call kim_set_model_numbering( & 777 model_driver_create_handle, KIM_NUMBERING_ONE_BASED, ierr) 778 if (ierr /= 0_c_int) then 779 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 780 "Unable to set numbering") 781 goto 42 782 end if 783 784 ! Store callback pointers in KIM object 785 call kim_set_routine_pointer( & 786 model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE, & 787 KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(compute), ierr) 788 call kim_set_routine_pointer( & 789 model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_CREATE, & 790 KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(compute_arguments_create), ierr2) 791 ierr = ierr + ierr2 792 call kim_set_routine_pointer( & 793 model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_DESTROY, & 794 KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(compute_arguments_destroy), ierr2) 795 ierr = ierr + ierr2 796 call kim_set_routine_pointer( & 797 model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_REFRESH, & 798 KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(refresh), ierr2) 799 ierr = ierr + ierr2 800 call kim_set_routine_pointer( & 801 model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_DESTROY, & 802 KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(destroy), ierr2) 803 ierr = ierr + ierr2 804 if (ierr /= 0_c_int) then 805 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 806 "Unable to store callback pointers") 807 goto 42 808 end if 809 810 ! Process parameter files 811 call kim_get_number_of_parameter_files( & 812 model_driver_create_handle, number_of_parameter_files) 813 if (number_of_parameter_files .ne. 1_c_int) then 814 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 815 "Wrong number of parameter files") 816 ierr = 1_c_int 817 goto 42 818 end if 819 820 ! Read in model parameters from parameter file 821 call kim_get_parameter_file_name( & 822 model_driver_create_handle, 1_c_int, parameter_file_name, ierr) 823 if (ierr /= 0_c_int) then 824 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 825 "Unable to get parameter file name") 826 ierr = 1_c_int 827 goto 42 828 end if 829 830 allocate(buf) 831 832 ! Read in model parameters from parameter file 833 ! 834 open(10,file=parameter_file_name,status="old") 835 read(10,*,iostat=ierr,err=100) in_species 836 read(10,*,iostat=ierr,err=100) in_cutoff 837 read(10,*,iostat=ierr,err=100) in_eps 838 read(10,*,iostat=ierr,err=100) in_sigma 839 close(10) 840 841 goto 200 842 100 continue 843 ! Reading parameters failed 844 ierr = 1_c_int 845 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 846 "Unable to read model parameters") 847 goto 42 848 849 200 continue 850 851 ! Convert parameters to requested units 852 call kim_convert_unit( & 853 KIM_LENGTH_UNIT_A, & 854 KIM_ENERGY_UNIT_EV, & 855 KIM_CHARGE_UNIT_E, & 856 KIM_TEMPERATURE_UNIT_K, & 857 KIM_TIME_UNIT_PS, & 858 requested_length_unit, & 859 requested_energy_unit, & 860 requested_charge_unit, & 861 requested_temperature_unit, & 862 requested_time_unit, & 863 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, convert_length, ierr) 864 if (ierr /= 0_c_int) then 865 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 866 "Failed to convert units") 867 goto 42 868 endif 869 870 in_cutoff = in_cutoff*convert_length 871 in_sigma = in_sigma*convert_length 872 873 call kim_convert_unit( & 874 KIM_LENGTH_UNIT_A, & 875 KIM_ENERGY_UNIT_EV, & 876 KIM_CHARGE_UNIT_E, & 877 KIM_TEMPERATURE_UNIT_K, & 878 KIM_TIME_UNIT_PS, & 879 requested_length_unit, & 880 requested_energy_unit, & 881 requested_charge_unit, & 882 requested_temperature_unit, & 883 requested_time_unit, & 884 0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, convert_energy, ierr) 885 if (ierr /= 0_c_int) then 886 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 887 "Failed to convert units") 888 goto 42 889 endif 890 891 in_eps = in_eps*convert_energy 892 893 ! Register units 894 call kim_set_units( & 895 model_driver_create_handle, & 896 requested_length_unit, & 897 requested_energy_unit, & 898 KIM_CHARGE_UNIT_UNUSED, & 899 KIM_TEMPERATURE_UNIT_UNUSED, & 900 KIM_TIME_UNIT_UNUSED, ierr) 901 if (ierr /= 0_c_int) then 902 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 903 "Unable to set units") 904 goto 42 905 end if 906 907 buf%influence_distance = in_cutoff 908 buf%cutoff(1) = in_cutoff 909 buf%cutsq(1) = in_cutoff*in_cutoff 910 buf%eps(1) = in_eps 911 buf%sigma(1) = in_sigma 912 buf%model_will_not_request_neighbors_of_noncontributing_particles(1) = 1_c_int 913 914 ! Compute energy at cutoff and store it in the buffer, as well 915 call calc_phi(in_eps, & 916 in_sigma, & 917 0.0_cd, & 918 in_cutoff, & 919 in_cutoff, energy_at_cutoff) 920 buf%shift = -energy_at_cutoff 921 922 ! Store model cutoff in KIM object 923 call kim_set_influence_distance_pointer( & 924 model_driver_create_handle, buf%influence_distance) 925 call kim_set_neighbor_list_pointers( & 926 model_driver_create_handle, 1_c_int, buf%cutoff, & 927 buf%model_will_not_request_neighbors_of_noncontributing_particles) 928 929 ! Register buffer in KIM API object 930 call kim_set_model_buffer_pointer( & 931 model_driver_create_handle, c_loc(buf)) 932 933 ! Register species 934 call kim_from_string(in_species, species_name) 935 936 call kim_set_species_code( & 937 model_driver_create_handle, species_name, speccode, ierr) 938 if (ierr /= 0_c_int) then 939 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 940 "Unable to set species code") 941 goto 42 942 end if 943 944 ! Set parameter pointers so they can be changed by the simulator 945 call kim_set_parameter_pointer( & 946 model_driver_create_handle, buf%cutoff, "cutoff", & 947 "Distance beyond which particles to not interact.", ierr) 948 if (ierr /= 0_c_int) then 949 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 950 "Failed to set parameter 'cutoff'.") 951 goto 42 952 endif 953 954 call kim_set_parameter_pointer( & 955 model_driver_create_handle, buf%eps, "epsilon", & 956 "Energy scaling coefficient (4*epsilon is the depth of the potential & 957 &well).", ierr) 958 if (ierr /= 0_c_int) then 959 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 960 "Failed to set parameter 'epsilon'.") 961 goto 42 962 endif 963 964 call kim_set_parameter_pointer( & 965 model_driver_create_handle, buf%sigma, "sigma", & 966 "Distance at which energy is equal to the 'shift' value added onto the & 967 &standard Lennard-Jones form.", ierr) 968 if (ierr /= 0_c_int) then 969 call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & 970 "Failed to set parameter 'sigma'.") 971 goto 42 972 endif 973 974 975 ierr = 0_c_int 976 42 continue 977 return 978 979end subroutine create 980