45 #include "KIM_API_status.h" 46 #define THIS_FILE_NAME __FILE__ 47 #define TRUEFALSE(TRUTH) merge(1,0,(TRUTH)) 51 use,
intrinsic :: iso_c_binding
61 integer(c_int),
parameter :: cd = c_double
62 integer(c_int),
parameter :: DIM = 3
63 integer(c_int),
parameter :: speccode = 1
66 real(c_double),
parameter :: model_cutsq =
model_cutoff**2
83 real(c_double),
parameter :: lj_epsilon = 0.0104_cd
84 real(c_double),
parameter :: lj_sigma = 3.40_cd
85 real(c_double),
parameter :: lj_cutnorm =
model_cutoff/lj_sigma
86 real(c_double),
parameter :: lj_a = 12.0_cd*lj_epsilon*(-26.0_cd &
87 + 7.0_cd*lj_cutnorm**6)/(lj_cutnorm**14 &
89 real(c_double),
parameter :: lj_b = 96.0_cd*lj_epsilon*(7.0_cd &
90 - 2.0_cd*lj_cutnorm**6)/(lj_cutnorm**13*lj_sigma)
91 real(c_double),
parameter :: lj_c = 28.0_cd*lj_epsilon*(-13.0_cd &
92 + 4.0_cd*lj_cutnorm**6)/(lj_cutnorm**12)
101 subroutine calc_phi(r,phi)
105 real(c_double),
intent(in) :: r
106 real(c_double),
intent(out) :: phi
109 real(c_double) rsq,sor,sor6,sor12
120 phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
123 end subroutine calc_phi
130 subroutine calc_phi_dphi(r,phi,dphi)
134 real(c_double),
intent(in) :: r
135 real(c_double),
intent(out) :: phi,dphi
138 real(c_double) rsq,sor,sor6,sor12
150 phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
151 dphi = 24.0_cd*lj_epsilon*(-2.0_cd*sor12+sor6)/r + 2.0_cd*lj_a*r + lj_b
154 end subroutine calc_phi_dphi
165 type(c_ptr),
intent(in) :: pkim
168 real(c_double) :: rij(dim)
169 real(c_double) :: r,rsqij,phi,dphi,deidr = 0.0_cd
170 integer(c_int) :: i,j,jj,numnei,part_ret,comp_force,comp_enepot,comp_virial, &
172 character (len=80) :: error_message
175 integer(c_int),
pointer :: n;
type(c_ptr) :: pn
176 real(c_double),
pointer :: energy;
type(c_ptr) :: penergy
177 real(c_double),
pointer :: coor(:,:);
type(c_ptr) :: pcoor
178 real(c_double),
pointer :: force(:,:);
type(c_ptr) :: pforce
179 real(c_double),
pointer :: enepot(:);
type(c_ptr) :: penepot
180 real(c_double),
pointer :: boxsidelengths(:);
type(c_ptr) :: pboxsidelengths
181 real(c_double),
pointer :: rij_list(:,:);
type(c_ptr) :: prij_list
182 integer(c_int),
pointer :: numcontrib;
type(c_ptr) :: pnumcontrib
183 integer(c_int),
pointer :: nei1part(:);
type(c_ptr) :: pnei1part
184 integer(c_int),
pointer :: particlespecies(:);
type(c_ptr) :: pparticlespecies
185 real(c_double),
pointer :: virial(:);
type(c_ptr) :: pvirial
186 character(len=KIM_KEY_STRING_LENGTH) nbc_method
187 integer(c_int) iterorloca
188 integer(c_int) halforfull
190 integer(c_int) numbercontrib
205 idum = kim_api_report_error(__line__, this_file_name, &
206 "kim_api_get_nbc_method", &
210 if (index(nbc_method,
"NEIGH_RVEC_H").eq.1)
then 213 elseif (index(nbc_method,
"NEIGH_RVEC_F").eq.1)
then 216 elseif (index(nbc_method,
"NEIGH_PURE_H").eq.1)
then 219 elseif (index(nbc_method,
"NEIGH_PURE_F").eq.1)
then 222 elseif (index(nbc_method,
"MI_OPBC_H").eq.1)
then 225 elseif (index(nbc_method,
"MI_OPBC_F").eq.1)
then 228 elseif (index(nbc_method,
"CLUSTER").eq.1)
then 233 idum = kim_api_report_error(__line__, this_file_name, &
247 idum = kim_api_report_error(__line__, this_file_name, &
248 "kim_api_get_neigh_mode", &
252 if (iterorloca.ne.1 .and. iterorloca.ne.2)
then 254 write(error_message,
'(a,i1)') &
255 'Unsupported IterOrLoca mode = ',iterorloca
256 idum = kim_api_report_error(__line__, this_file_name, &
268 "energy", comp_energy, 1, &
269 "forces", comp_force, 1, &
270 "particleEnergy", comp_enepot, 1, &
271 "virial", comp_virial, 1)
273 idum = kim_api_report_error(__line__, this_file_name, &
281 "numberOfParticles", pn, 1, &
282 "particleSpecies", pparticlespecies,1, &
283 "coordinates", pcoor, 1, &
284 "numberContributingParticles", pnumcontrib, truefalse(halforfull.eq.1), &
285 "boxSideLengths", pboxsidelengths, truefalse(nbc.eq.2), &
286 "energy", penergy, truefalse(comp_energy.eq.1), &
287 "forces", pforce, truefalse(comp_force.eq.1), &
288 "particleEnergy", penepot, truefalse(comp_enepot.eq.1), &
289 "virial", pvirial, truefalse(comp_virial.eq.1))
291 idum = kim_api_report_error(__line__, this_file_name, &
296 call c_f_pointer(pn, n)
297 call c_f_pointer(pparticlespecies, particlespecies, [n])
298 call c_f_pointer(pcoor, coor, [dim,n])
299 if (halforfull.eq.1)
call c_f_pointer(pnumcontrib, numcontrib)
300 if (nbc.eq.2)
call c_f_pointer(pboxsidelengths, boxsidelengths, [dim])
301 if (comp_energy.eq.1)
call c_f_pointer(penergy, energy)
302 if (comp_force.eq.1)
call c_f_pointer(pforce, force, [dim,n])
303 if (comp_enepot.eq.1)
call c_f_pointer(penepot, enepot, [n])
304 if (comp_virial.eq.1)
call c_f_pointer(pvirial, virial, [6])
306 if (halforfull.eq.1)
then 308 numbercontrib = numcontrib
318 if (particlespecies(i).ne.speccode)
then 319 idum = kim_api_report_error(__line__, this_file_name, &
320 "Unexpected species detected", &
329 if (comp_enepot.eq.1) enepot = 0.0_cd
330 if (comp_energy.eq.1) energy = 0.0_cd
331 if (comp_force.eq.1) force = 0.0_cd
332 if (comp_virial.eq.1) virial = 0.0_cd
337 allocate( nei1part(n) )
342 if (iterorloca.eq.1)
then 347 idum = kim_api_report_error(__line__, this_file_name, &
365 if (iterorloca.eq.1)
then 373 idum = kim_api_report_error(__line__, this_file_name, &
374 "kim_api_get_neigh", &
387 nei1part(1:numnei) = (/ (i+jj, jj = 1,numnei) /)
394 idum = kim_api_report_error(__line__, this_file_name, &
395 "kim_api_get_neigh", &
403 if (nbc.ne.3)
call c_f_pointer(pnei1part, nei1part, [numnei])
404 if (nbc.eq.0)
call c_f_pointer(prij_list, rij_list, [dim,numnei])
415 rij(:) = coor(:,j) - coor(:,i)
417 rij(:) = rij_list(:,jj)
423 where ( abs(rij) .gt. 0.5_cd*boxsidelengths )
425 rij = rij - sign(boxsidelengths,rij)
431 rsqij = dot_product(rij,rij)
432 if ( rsqij .lt. model_cutsq )
then 435 if (comp_force.eq.1.or.comp_virial.eq.1)
then 436 call calc_phi_dphi(r,phi,dphi)
438 if ((halforfull.eq.1) .and. &
439 (j .le. numbercontrib))
then 450 if (comp_enepot.eq.1)
then 451 enepot(i) = enepot(i) + 0.5_cd*phi
452 if ((halforfull.eq.1) .and. &
453 (j .le. numbercontrib)) &
454 enepot(j) = enepot(j) + 0.5_cd*phi
456 if (comp_energy.eq.1)
then 457 if ((halforfull.eq.1) .and. &
458 (j .le. numbercontrib))
then 459 energy = energy + phi
461 energy = energy + 0.5_cd*phi
467 if (comp_virial.eq.1)
then 468 virial(1) = virial(1) + rij(1)*rij(1)*deidr/r
469 virial(2) = virial(2) + rij(2)*rij(2)*deidr/r
470 virial(3) = virial(3) + rij(3)*rij(3)*deidr/r
471 virial(4) = virial(4) + rij(2)*rij(3)*deidr/r
472 virial(5) = virial(5) + rij(1)*rij(3)*deidr/r
473 virial(6) = virial(6) + rij(1)*rij(2)*deidr/r
478 if (comp_force.eq.1)
then 479 force(:,i) = force(:,i) + deidr*rij/r
480 force(:,j) = force(:,j) - deidr*rij/r
491 if (nbc.eq.3)
deallocate( nei1part )
507 integer(c_int) function model_init(pkim) bind(c)
508 use,
intrinsic :: iso_c_binding
514 type(c_ptr),
intent(in) :: pkim
517 integer(c_int),
parameter :: one=1
518 integer(c_int) ier, idum
521 real(c_double),
pointer :: cutoff;
type(c_ptr) :: pcutoff
525 if (ier.lt.kim_status_ok)
then 526 idum = kim_api_report_error(__line__, this_file_name, &
527 "kim_api_set_method", ier)
532 pcutoff = kim_api_get_data(pkim,
"cutoff",ier)
533 if (ier.lt.kim_status_ok)
then 534 idum = kim_api_report_error(__line__, this_file_name, &
535 "kim_api_get_data", ier)
538 call c_f_pointer(pcutoff, cutoff)
546 end function model_init
real(c_double), parameter, public model_cutoff
integer(c_int) function, public compute_energy_forces(pkim)