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