KIM API V2
ex_model_Ar_P_MLJ_F03.F03
Go to the documentation of this file.
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--2018, 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 ex_model_Ar_P_MLJ_F03
35 !**
36 !** Modified Lennard-Jones pair potential (with smooth cutoff) model for Ar
37 !**
38 !** Reference: Ashcroft and Mermin
39 !**
40 !** Language: Fortran 2003
41 !**
42 !****************************************************************************
43 
44 
46 
47 use, intrinsic :: iso_c_binding
49 implicit none
50 
51 save
52 private
53 public &!Compute_Energy_Forces, &
54  model_cutoff, &
55  speccode, &
56  buffer_type
57 
58 ! Below are the definitions and values of all Model parameters
59 integer(c_int), parameter :: cd = c_double ! used for literal constants
60 integer(c_int), parameter :: DIM = 3 ! dimensionality of space
61 integer(c_int), parameter :: speccode = 1 ! internal species code
62 real(c_double), parameter :: model_cutoff = 8.15_cd ! cutoff radius
63  ! in angstroms
64 real(c_double), parameter :: model_cutsq = model_cutoff**2
65 
66 !-------------------------------------------------------------------------------
67 ! Below are the definitions and values of all additional model parameters
68 !
69 ! Recall that the Fortran 2003 format for declaring parameters is as follows:
70 !
71 ! integer(c_int), parameter :: parname = value ! This defines an integer
72 ! ! parameter called `parname'
73 ! ! with a value equal to
74 ! ! `value' (a number)
75 !
76 ! real(c_double), parameter :: parname = value ! This defines a real(c_double)
77 ! ! parameter called `parname'
78 ! ! with a value equal to
79 ! ! `value' (a number)
80 !-------------------------------------------------------------------------------
81 real(c_double), parameter :: lj_epsilon = 0.0104_cd
82 real(c_double), parameter :: lj_sigma = 3.40_cd
83 real(c_double), parameter :: lj_cutnorm = model_cutoff/lj_sigma
84 real(c_double), parameter :: lj_A = 12.0_cd*lj_epsilon*(-26.0_cd &
85  + 7.0_cd*lj_cutnorm**6)/(lj_cutnorm**14 &
86  *lj_sigma**2)
87 real(c_double), parameter :: lj_B = 96.0_cd*lj_epsilon*(7.0_cd &
88  - 2.0_cd*lj_cutnorm**6)/(lj_cutnorm**13*lj_sigma)
89 real(c_double), parameter :: lj_C = 28.0_cd*lj_epsilon*(-13.0_cd &
90  + 4.0_cd*lj_cutnorm**6)/(lj_cutnorm**12)
91 
92 type, bind(c) :: buffer_type
93  real(c_double) :: influence_distance
94  real(c_double) :: cutoff(1)
95 end type buffer_type
96 
97 contains
98 
99 !-------------------------------------------------------------------------------
100 !
101 ! Calculate pair potential phi(r)
102 !
103 !-------------------------------------------------------------------------------
104 subroutine calc_phi(r,phi)
105 implicit none
106 
107 !-- Transferred variables
108 real(c_double), intent(in) :: r
109 real(c_double), intent(out) :: phi
110 
111 !-- Local variables
112 real(c_double) rsq,sor,sor6,sor12
113 
114 rsq = r*r ! r^2
115 sor = lj_sigma/r ! (sig/r)
116 sor6 = sor*sor*sor !
117 sor6 = sor6*sor6 ! (sig/r)^6
118 sor12= sor6*sor6 ! (sig/r)^12
119 if (r .gt. model_cutoff) then
120  ! Argument exceeds cutoff radius
121  phi = 0.0_cd
122 else
123  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
124 endif
125 
126 end subroutine calc_phi
127 
128 !-------------------------------------------------------------------------------
129 !
130 ! Calculate pair potential phi(r) and its derivative dphi(r)
131 !
132 !-------------------------------------------------------------------------------
133 subroutine calc_phi_dphi(r,phi,dphi)
134 implicit none
135 
136 !-- Transferred variables
137 real(c_double), intent(in) :: r
138 real(c_double), intent(out) :: phi,dphi
139 
140 !-- Local variables
141 real(c_double) rsq,sor,sor6,sor12
142 
143 rsq = r*r ! r^2
144 sor = lj_sigma/r ! (sig/r)
145 sor6 = sor*sor*sor !
146 sor6 = sor6*sor6 ! (sig/r)^6
147 sor12= sor6*sor6 ! (sig/r)^12
148 if (r .gt. model_cutoff) then
149  ! Argument exceeds cutoff radius
150  phi = 0.0_cd
151  dphi = 0.0_cd
152 else
153  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
154  dphi = 24.0_cd*lj_epsilon*(-2.0_cd*sor12+sor6)/r + 2.0_cd*lj_a*r + lj_b
155 endif
156 
157 end subroutine calc_phi_dphi
158 
159 !!-------------------------------------------------------------------------------
160 !!
161 !! Compute energy and forces on particles from the positions.
162 !!
163 !!-------------------------------------------------------------------------------
164 !integer(c_int) function Compute_Energy_Forces(pkim) bind(c)
165 !implicit none
166 !
167 !!-- Transferred variables
168 !type(c_ptr), intent(in) :: pkim
169 !
170 !!-- Local variables
171 !real(c_double) :: Rij(DIM)
172 !real(c_double) :: r,Rsqij,phi,dphi,dEidr = 0.0_cd
173 !integer(c_int) :: i,j,jj,numnei,part_ret,comp_force,comp_enepot,comp_virial, &
174 ! comp_energy
175 !character (len=80) :: error_message
176 !
177 !!-- KIM variables
178 !integer(c_int), pointer :: N; type(c_ptr) :: pN
179 !real(c_double), pointer :: energy; type(c_ptr) :: penergy
180 !real(c_double), pointer :: coor(:,:); type(c_ptr) :: pcoor
181 !real(c_double), pointer :: force(:,:); type(c_ptr) :: pforce
182 !real(c_double), pointer :: enepot(:); type(c_ptr) :: penepot
183 !real(c_double), pointer :: Rij_list(:,:); type(c_ptr) :: pRij_list
184 !integer(c_int), pointer :: nei1part(:); type(c_ptr) :: pnei1part
185 !integer(c_int), pointer :: particleSpecies(:);type(c_ptr) :: pparticleSpecies
186 !real(c_double), pointer :: virial(:); type(c_ptr) :: pvirial
187 !integer(c_int) idum
188 !
189 !
190 !! Check to see if we have been asked to compute the forces, energyperpart,
191 !! energy and virial
192 !!
193 !call kim_api_getm_compute(pkim, Compute_Energy_Forces, &
194 ! "energy", comp_energy, 1, &
195 ! "forces", comp_force, 1, &
196 ! "particleEnergy", comp_enepot, 1, &
197 ! "virial", comp_virial, 1)
198 !if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
199 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
200 ! "kim_api_getm_compute", Compute_Energy_Forces)
201 ! return
202 !endif
203 !
204 !! Unpack data from KIM object
205 !!
206 !call kim_api_getm_data(pkim, Compute_Energy_Forces, &
207 ! "numberOfParticles", pN, 1, &
208 ! "particleSpecies", pparticleSpecies,1, &
209 ! "coordinates", pcoor, 1, &
210 ! "energy", penergy, TRUEFALSE(comp_energy.eq.1), &
211 ! "forces", pforce, TRUEFALSE(comp_force.eq.1), &
212 ! "particleEnergy", penepot, TRUEFALSE(comp_enepot.eq.1), &
213 ! "virial", pvirial, TRUEFALSE(comp_virial.eq.1))
214 !if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
215 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
216 ! "kim_api_getm_data_f", Compute_Energy_Forces)
217 ! return
218 !endif
219 !
220 !call c_f_pointer(pN, N)
221 !call c_f_pointer(pparticleSpecies, particleSpecies, [N])
222 !call c_f_pointer(pcoor, coor, [DIM,N])
223 !if (comp_energy.eq.1) call c_f_pointer(penergy, energy)
224 !if (comp_force.eq.1) call c_f_pointer(pforce, force, [DIM,N])
225 !if (comp_enepot.eq.1) call c_f_pointer(penepot, enepot, [N])
226 !if (comp_virial.eq.1) call c_f_pointer(pvirial, virial, [6])
227 !
228 !
229 !! Check to be sure that the species are correct
230 !!
231 !Compute_Energy_Forces = KIM_STATUS_FAIL ! assume an error
232 !do i = 1,N
233 ! if (particleSpecies(i).ne.speccode) then
234 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
235 ! "Unexpected species detected", &
236 ! Compute_Energy_Forces)
237 ! return
238 ! endif
239 !enddo
240 !Compute_Energy_Forces = KIM_STATUS_OK ! everything is ok
241 !
242 !! Initialize potential energies, forces, virial term
243 !!
244 !if (comp_enepot.eq.1) enepot = 0.0_cd
245 !if (comp_energy.eq.1) energy = 0.0_cd
246 !if (comp_force.eq.1) force = 0.0_cd
247 !if (comp_virial.eq.1) virial = 0.0_cd
248 !
249 !
250 !!
251 !! Compute energy and forces
252 !!
253 !
254 !! Loop over particles and compute energy and forces
255 !!
256 !do i=1,N
257 ! Compute_Energy_Forces = kim_api_get_neigh(pkim,1,i,part_ret,numnei, &
258 ! pnei1part,pRij_list)
259 ! if (Compute_Energy_Forces.ne.KIM_STATUS_OK) then
260 ! ! some sort of problem, exit
261 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
262 ! "kim_api_get_neigh", &
263 ! Compute_Energy_Forces)
264 ! Compute_Energy_Forces = KIM_STATUS_FAIL
265 ! return
266 ! endif
267 !
268 ! call c_f_pointer(pnei1part, nei1part, [numnei])
269 !
270 ! ! Loop over the neighbors of particle i
271 ! !
272 ! do jj = 1, numnei
273 !
274 ! j = nei1part(jj) ! get neighbor ID
275 !
276 ! ! compute relative position vector
277 ! !
278 ! Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
279 !
280 ! ! compute energy and forces
281 ! !
282 ! Rsqij = dot_product(Rij,Rij) ! compute square distance
283 ! if ( Rsqij .lt. model_cutsq ) then ! particles are interacting?
284 !
285 ! r = sqrt(Rsqij) ! compute distance
286 ! if (comp_force.eq.1.or.comp_virial.eq.1) then
287 ! call calc_phi_dphi(r,phi,dphi) ! compute pair potential
288 ! ! and it derivative
289 ! dEidr = 0.5_cd*dphi
290 ! else
291 ! call calc_phi(r,phi) ! compute just pair potential
292 ! endif
293 !
294 ! ! contribution to energy
295 ! !
296 ! if (comp_enepot.eq.1) then
297 ! enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
298 ! endif
299 ! if (comp_energy.eq.1) then
300 ! energy = energy + 0.5_cd*phi
301 ! endif
302 !
303 ! ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
304 ! !
305 ! if (comp_virial.eq.1) then
306 ! virial(1) = virial(1) + Rij(1)*Rij(1)*dEidr/r
307 ! virial(2) = virial(2) + Rij(2)*Rij(2)*dEidr/r
308 ! virial(3) = virial(3) + Rij(3)*Rij(3)*dEidr/r
309 ! virial(4) = virial(4) + Rij(2)*Rij(3)*dEidr/r
310 ! virial(5) = virial(5) + Rij(1)*Rij(3)*dEidr/r
311 ! virial(6) = virial(6) + Rij(1)*Rij(2)*dEidr/r
312 ! endif
313 !
314 ! ! contribution to forces
315 ! !
316 ! if (comp_force.eq.1) then
317 ! force(:,i) = force(:,i) + dEidr*Rij/r ! accumulate force on i
318 ! force(:,j) = force(:,j) - dEidr*Rij/r ! accumulate force on j
319 ! endif
320 !
321 ! endif
322 !
323 ! enddo ! loop on jj
324 !
325 !enddo
326 !
327 !! Everything is great
328 !!
329 !Compute_Energy_Forces = KIM_STATUS_OK
330 !return
331 !
332 !end function Compute_Energy_Forces
333 
334 end module ex_model_ar_p_mlj_f03
335 
336 !-------------------------------------------------------------------------------
337 !
338 ! Model create routine (REQUIRED)
339 !
340 !-------------------------------------------------------------------------------
341 #include "kim_model_create_log_macros.fd"
342 subroutine model_create_routine(model_create_handle, requested_length_unit, &
343  requested_energy_unit, requested_charge_unit, requested_temperature_unit, &
344  requested_time_unit, ierr) bind(c)
345 use, intrinsic :: iso_c_binding
348 implicit none
349 
350 !-- Transferred variables
351 type(kim_model_create_handle_type), intent(inout) :: model_create_handle
352 type(kim_length_unit_type), intent(in) :: requested_length_unit
353 type(kim_energy_unit_type), intent(in) :: requested_energy_unit
354 type(kim_charge_unit_type), intent(in) :: requested_charge_unit
355 type(kim_temperature_unit_type), intent(in) :: requested_temperature_unit
356 type(kim_time_unit_type), intent(in) :: requested_time_unit
357 integer(c_int), intent(out) :: ierr
358 
359 !-- KIM variables
360 integer(c_int) :: ierr2
361 type(buffer_type), pointer :: buf
362 character(kind=c_char, len=100000) model_create_string
363 
364 kim_log_file = __file__
365 
366 ierr = 0
367 ierr2 = 0
368 
369 ! set units
370 call kim_model_create_set_units(model_create_handle, &
371  kim_length_unit_a, &
372  kim_energy_unit_ev, &
373  kim_charge_unit_unused, &
374  kim_temperature_unit_unused, &
375  kim_time_unit_unused, &
376  ierr2)
377 ierr = ierr + ierr2
378 
379 ! register species
380 call kim_model_create_set_species_code(model_create_handle, &
381  kim_species_name_ar, speccode, ierr2)
382 ierr = ierr + ierr2
383 
384 ! register numbering
385 call kim_model_create_set_model_numbering(model_create_handle, &
386  kim_numbering_one_based, ierr2);
387 ierr = ierr + ierr2
388 
389 ! register function pointers
390 call kim_model_create_set_compute_pointer(model_create_handle, &
391  kim_language_name_fortran, c_funloc(kim_model_create_string), ierr2)
392 ierr = ierr + ierr2
393 call kim_model_create_set_compute_arguments_create_pointer( &
394  model_create_handle, kim_language_name_fortran, &
395  c_funloc(kim_model_create_string), ierr2)
396 ierr = ierr + ierr2
397 call kim_model_create_set_compute_arguments_destroy_pointer( &
398  model_create_handle, kim_language_name_fortran, &
399  c_funloc(kim_model_create_string), ierr2)
400 ierr = ierr + ierr2
401 call kim_model_create_set_destroy_pointer(model_create_handle, &
402  kim_language_name_fortran, c_funloc(kim_model_create_string), ierr2)
403 ierr = ierr + ierr2
404 call kim_model_create_set_refresh_pointer( &
405  model_create_handle, kim_language_name_fortran, &
406  c_funloc(kim_model_create_string), ierr2)
407 ierr = ierr + ierr2
408 
409 ! allocate buffer
410 allocate( buf )
411 
412 ! store model buffer in KIM object
413 call kim_model_create_set_model_buffer_pointer(model_create_handle, &
414  c_loc(buf))
415 
416 ! set buffer values
417 buf%influence_distance = model_cutoff
418 buf%cutoff = model_cutoff
419 
420 ! register influence distance
421 call kim_model_create_set_influence_distance_pointer( &
422  model_create_handle, buf%influence_distance)
423 
424 ! register cutoff
425 call kim_model_create_set_neighbor_list_cutoffs_pointer(model_create_handle, &
426  1, buf%cutoff)
427 
428 if (ierr /= 0) then
429  ierr = 1
430  deallocate( buf )
431  kim_log_message = "Unable to successfully initialize model"
432  log_error()
433 endif
434 
435 call kim_model_create_string(model_create_handle, model_create_string)
436 print '(A)', trim(model_create_string)
437 
438 return
439 
440 end subroutine model_create_routine