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
64 real(c_double),
parameter :: model_cutoff = 8.15_cd
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
116 if (r .gt. model_cutoff)
then 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
145 if (r .gt. model_cutoff)
then 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 :: Rij_list(:,:);
type(c_ptr) :: pRij_list
181 integer(c_int),
pointer :: nei1part(:);
type(c_ptr) :: pnei1part
182 integer(c_int),
pointer :: particleSpecies(:);
type(c_ptr) :: pparticleSpecies
183 real(c_double),
pointer :: virial(:);
type(c_ptr) :: pvirial
191 "energy", comp_energy, 1, &
192 "forces", comp_force, 1, &
193 "particleEnergy", comp_enepot, 1, &
194 "virial", comp_virial, 1)
196 idum = kim_api_report_error(__line__, this_file_name, &
204 "numberOfParticles", pn, 1, &
205 "particleSpecies", pparticlespecies,1, &
206 "coordinates", pcoor, 1, &
207 "energy", penergy, truefalse(comp_energy.eq.1), &
208 "forces", pforce, truefalse(comp_force.eq.1), &
209 "particleEnergy", penepot, truefalse(comp_enepot.eq.1), &
210 "virial", pvirial, truefalse(comp_virial.eq.1))
212 idum = kim_api_report_error(__line__, this_file_name, &
217 call c_f_pointer(pn, n)
218 call c_f_pointer(pparticlespecies, particlespecies, [n])
219 call c_f_pointer(pcoor, coor, [dim,n])
220 if (comp_energy.eq.1)
call c_f_pointer(penergy, energy)
221 if (comp_force.eq.1)
call c_f_pointer(pforce, force, [dim,n])
222 if (comp_enepot.eq.1)
call c_f_pointer(penepot, enepot, [n])
223 if (comp_virial.eq.1)
call c_f_pointer(pvirial, virial, [6])
230 if (particlespecies(i).ne.speccode)
then 231 idum = kim_api_report_error(__line__, this_file_name, &
232 "Unexpected species detected", &
241 if (comp_enepot.eq.1) enepot = 0.0_cd
242 if (comp_energy.eq.1) energy = 0.0_cd
243 if (comp_force.eq.1) force = 0.0_cd
244 if (comp_virial.eq.1) virial = 0.0_cd
258 idum = kim_api_report_error(__line__, this_file_name, &
259 "kim_api_get_neigh", &
265 call c_f_pointer(pnei1part, nei1part, [numnei])
275 rij(:) = coor(:,j) - coor(:,i)
279 rsqij = dot_product(rij,rij)
280 if ( rsqij .lt. model_cutsq )
then 283 if (comp_force.eq.1.or.comp_virial.eq.1)
then 284 call calc_phi_dphi(r,phi,dphi)
293 if (comp_enepot.eq.1)
then 294 enepot(i) = enepot(i) + 0.5_cd*phi
296 if (comp_energy.eq.1)
then 297 energy = energy + 0.5_cd*phi
302 if (comp_virial.eq.1)
then 303 virial(1) = virial(1) + rij(1)*rij(1)*deidr/r
304 virial(2) = virial(2) + rij(2)*rij(2)*deidr/r
305 virial(3) = virial(3) + rij(3)*rij(3)*deidr/r
306 virial(4) = virial(4) + rij(2)*rij(3)*deidr/r
307 virial(5) = virial(5) + rij(1)*rij(3)*deidr/r
308 virial(6) = virial(6) + rij(1)*rij(2)*deidr/r
313 if (comp_force.eq.1)
then 314 force(:,i) = force(:,i) + deidr*rij/r
315 force(:,j) = force(:,j) - deidr*rij/r
338 integer(c_int) function model_init(pkim) bind(c)
339 use,
intrinsic :: iso_c_binding
345 type(c_ptr),
intent(in) :: pkim
348 integer(c_int),
parameter :: one=1
349 integer(c_int) ier, idum
352 real(c_double),
pointer :: cutoff;
type(c_ptr) :: pcutoff
356 if (ier.lt.kim_status_ok)
then 357 idum = kim_api_report_error(__line__, this_file_name, &
358 "kim_api_set_method", ier)
363 pcutoff = kim_api_get_data(pkim,
"cutoff",ier)
364 if (ier.lt.kim_status_ok)
then 365 idum = kim_api_report_error(__line__, this_file_name, &
366 "kim_api_get_data", ier)
369 call c_f_pointer(pcutoff, cutoff)
370 cutoff = model_cutoff
377 end function model_init
real(c_double), parameter, public model_cutoff
integer(c_int) function, public compute_energy_forces(pkim)