ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/DarkSide/eam.F90
Revision: 578
Committed: Fri Aug 26 21:30:41 2005 UTC (19 years, 11 months ago) by chrisfen
File size: 32610 byte(s)
Log Message:
added some probably nonfunctional get*cut routines

File Contents

# Content
1 !!
2 !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 !!
4 !! The University of Notre Dame grants you ("Licensee") a
5 !! non-exclusive, royalty free, license to use, modify and
6 !! redistribute this software in source and binary code form, provided
7 !! that the following conditions are met:
8 !!
9 !! 1. Acknowledgement of the program authors must be made in any
10 !! publication of scientific results based in part on use of the
11 !! program. An acceptable form of acknowledgement is citation of
12 !! the article in which the program was described (Matthew
13 !! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 !! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 !! Parallel Simulation Engine for Molecular Dynamics,"
16 !! J. Comput. Chem. 26, pp. 252-271 (2005))
17 !!
18 !! 2. Redistributions of source code must retain the above copyright
19 !! notice, this list of conditions and the following disclaimer.
20 !!
21 !! 3. Redistributions in binary form must reproduce the above copyright
22 !! notice, this list of conditions and the following disclaimer in the
23 !! documentation and/or other materials provided with the
24 !! distribution.
25 !!
26 !! This software is provided "AS IS," without a warranty of any
27 !! kind. All express or implied conditions, representations and
28 !! warranties, including any implied warranty of merchantability,
29 !! fitness for a particular purpose or non-infringement, are hereby
30 !! excluded. The University of Notre Dame and its licensors shall not
31 !! be liable for any damages suffered by licensee as a result of
32 !! using, modifying or distributing the software or its
33 !! derivatives. In no event will the University of Notre Dame or its
34 !! licensors be liable for any lost revenue, profit or data, or for
35 !! direct, indirect, special, consequential, incidental or punitive
36 !! damages, however caused and regardless of the theory of liability,
37 !! arising out of the use of or inability to use software, even if the
38 !! University of Notre Dame has been advised of the possibility of
39 !! such damages.
40 !!
41
42 module eam
43 use simulation
44 use force_globals
45 use status
46 use atype_module
47 use vector_class
48 #ifdef IS_MPI
49 use mpiSimulation
50 #endif
51 implicit none
52 PRIVATE
53
54 INTEGER, PARAMETER :: DP = selected_real_kind(15)
55
56 logical, save :: EAM_FF_initialized = .false.
57 integer, save :: EAM_Mixing_Policy
58 real(kind = dp), save :: EAM_rcut
59 logical, save :: haveRcut = .false.
60
61 character(len = statusMsgSize) :: errMesg
62 integer :: eam_err
63
64 character(len = 200) :: errMsg
65 character(len=*), parameter :: RoutineName = "EAM MODULE"
66 !! Logical that determines if eam arrays should be zeroed
67 logical :: cleanme = .true.
68 logical :: nmflag = .false.
69
70
71 type, private :: EAMtype
72 integer :: eam_atype
73 real( kind = DP ) :: eam_dr
74 integer :: eam_nr
75 integer :: eam_nrho
76 real( kind = DP ) :: eam_lattice
77 real( kind = DP ) :: eam_drho
78 real( kind = DP ) :: eam_rcut
79 integer :: eam_atype_map
80
81 real( kind = DP ), pointer, dimension(:) :: eam_rvals => null()
82 real( kind = DP ), pointer, dimension(:) :: eam_rhovals => null()
83 real( kind = DP ), pointer, dimension(:) :: eam_F_rho => null()
84 real( kind = DP ), pointer, dimension(:) :: eam_Z_r => null()
85 real( kind = DP ), pointer, dimension(:) :: eam_rho_r => null()
86 real( kind = DP ), pointer, dimension(:) :: eam_phi_r => null()
87 real( kind = DP ), pointer, dimension(:) :: eam_F_rho_pp => null()
88 real( kind = DP ), pointer, dimension(:) :: eam_Z_r_pp => null()
89 real( kind = DP ), pointer, dimension(:) :: eam_rho_r_pp => null()
90 real( kind = DP ), pointer, dimension(:) :: eam_phi_r_pp => null()
91 end type EAMtype
92
93
94 !! Arrays for derivatives used in force calculation
95 real( kind = dp), dimension(:), allocatable :: frho
96 real( kind = dp), dimension(:), allocatable :: rho
97
98 real( kind = dp), dimension(:), allocatable :: dfrhodrho
99 real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho
100
101
102 !! Arrays for MPI storage
103 #ifdef IS_MPI
104 real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col
105 real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row
106 real( kind = dp),save, dimension(:), allocatable :: frho_row
107 real( kind = dp),save, dimension(:), allocatable :: frho_col
108 real( kind = dp),save, dimension(:), allocatable :: rho_row
109 real( kind = dp),save, dimension(:), allocatable :: rho_col
110 real( kind = dp),save, dimension(:), allocatable :: rho_tmp
111 real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col
112 real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row
113 #endif
114
115 type, private :: EAMTypeList
116 integer :: n_eam_types = 0
117 integer :: currentAddition = 0
118
119 type (EAMtype), pointer :: EAMParams(:) => null()
120 integer, pointer :: atidToEAMType(:) => null()
121 end type EAMTypeList
122
123
124 type (eamTypeList), save :: EAMList
125
126 !! standard eam stuff
127
128
129 public :: init_EAM_FF
130 public :: setCutoffEAM
131 public :: do_eam_pair
132 public :: newEAMtype
133 public :: calc_eam_prepair_rho
134 public :: calc_eam_preforce_Frho
135 public :: clean_EAM
136 public :: destroyEAMTypes
137 public :: getEAMCut
138
139 contains
140
141
142 subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
143 eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
144 eam_ident,status)
145 real (kind = dp ) :: lattice_constant
146 integer :: eam_nrho
147 real (kind = dp ) :: eam_drho
148 integer :: eam_nr
149 real (kind = dp ) :: eam_dr
150 real (kind = dp ) :: rcut
151 real (kind = dp ), dimension(eam_nr) :: eam_Z_r
152 real (kind = dp ), dimension(eam_nr) :: eam_rho_r
153 real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
154 integer :: eam_ident
155 integer :: status
156
157 integer :: nAtypes,nEAMTypes,myATID
158 integer :: maxVals
159 integer :: alloc_stat
160 integer :: current
161 integer,pointer :: Matchlist(:) => null()
162
163 status = 0
164
165
166 !! Assume that atypes has already been set and get the total number of types in atypes
167 !! Also assume that every member of atypes is a EAM model.
168
169
170 ! check to see if this is the first time into
171 if (.not.associated(EAMList%EAMParams)) then
172 call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
173 EAMList%n_eam_types = nEAMtypes
174 allocate(EAMList%EAMParams(nEAMTypes))
175 nAtypes = getSize(atypes)
176 allocate(EAMList%atidToEAMType(nAtypes))
177 end if
178
179 EAMList%currentAddition = EAMList%currentAddition + 1
180 current = EAMList%currentAddition
181
182 myATID = getFirstMatchingElement(atypes, "c_ident", eam_ident)
183 EAMList%atidToEAMType(myATID) = current
184
185 call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
186 if (alloc_stat /= 0) then
187 status = -1
188 return
189 end if
190
191
192 EAMList%EAMParams(current)%eam_atype = eam_ident
193 EAMList%EAMParams(current)%eam_lattice = lattice_constant
194 EAMList%EAMParams(current)%eam_nrho = eam_nrho
195 EAMList%EAMParams(current)%eam_drho = eam_drho
196 EAMList%EAMParams(current)%eam_nr = eam_nr
197 EAMList%EAMParams(current)%eam_dr = eam_dr
198 EAMList%EAMParams(current)%eam_rcut = rcut
199 EAMList%EAMParams(current)%eam_Z_r = eam_Z_r
200 EAMList%EAMParams(current)%eam_rho_r = eam_rho_r
201 EAMList%EAMParams(current)%eam_F_rho = eam_F_rho
202
203 end subroutine newEAMtype
204
205
206 ! kills all eam types entered and sets simulation to uninitalized
207 subroutine destroyEAMtypes()
208 integer :: i
209 type(EAMType), pointer :: tempEAMType=>null()
210
211 do i = 1, EAMList%n_eam_types
212 tempEAMType => eamList%EAMParams(i)
213 call deallocate_EAMType(tempEAMType)
214 end do
215 if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
216 eamList%EAMParams => null()
217
218 eamList%n_eam_types = 0
219 eamList%currentAddition = 0
220
221 end subroutine destroyEAMtypes
222
223 function getEAMCut(atomID) result(cutValue)
224 integer, intent(in) :: atomID
225 real(kind=dp) :: cutValue
226
227 cutValue = eamList%EAMParams(atomID)%eam_rcut
228
229 end function getEAMCut
230
231 subroutine init_EAM_FF(status)
232 integer :: status
233 integer :: i,j
234 real(kind=dp) :: current_rcut_max
235 integer :: alloc_stat
236 integer :: number_r, number_rho
237
238
239 status = 0
240 if (EAMList%currentAddition == 0) then
241 call handleError("init_EAM_FF","No members in EAMList")
242 status = -1
243 return
244 end if
245
246
247 do i = 1, EAMList%currentAddition
248
249 ! Build array of r values
250
251 do j = 1,EAMList%EAMParams(i)%eam_nr
252 EAMList%EAMParams(i)%eam_rvals(j) = &
253 real(j-1,kind=dp)* &
254 EAMList%EAMParams(i)%eam_dr
255 end do
256 ! Build array of rho values
257 do j = 1,EAMList%EAMParams(i)%eam_nrho
258 EAMList%EAMParams(i)%eam_rhovals(j) = &
259 real(j-1,kind=dp)* &
260 EAMList%EAMParams(i)%eam_drho
261 end do
262 ! convert from eV to kcal / mol:
263 EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP
264
265 ! precompute the pair potential and get it into kcal / mol:
266 EAMList%EAMParams(i)%eam_phi_r(1) = 0.0E0_DP
267 do j = 2, EAMList%EAMParams(i)%eam_nr
268 EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
269 EAMList%EAMParams(i)%eam_phi_r(j) = EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
270 enddo
271 end do
272
273
274 do i = 1, EAMList%currentAddition
275 number_r = EAMList%EAMParams(i)%eam_nr
276 number_rho = EAMList%EAMParams(i)%eam_nrho
277
278 call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
279 EAMList%EAMParams(i)%eam_rho_r, &
280 EAMList%EAMParams(i)%eam_rho_r_pp, &
281 0.0E0_DP, 0.0E0_DP, 'N')
282 call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
283 EAMList%EAMParams(i)%eam_Z_r, &
284 EAMList%EAMParams(i)%eam_Z_r_pp, &
285 0.0E0_DP, 0.0E0_DP, 'N')
286 call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
287 EAMList%EAMParams(i)%eam_F_rho, &
288 EAMList%EAMParams(i)%eam_F_rho_pp, &
289 0.0E0_DP, 0.0E0_DP, 'N')
290 call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
291 EAMList%EAMParams(i)%eam_phi_r, &
292 EAMList%EAMParams(i)%eam_phi_r_pp, &
293 0.0E0_DP, 0.0E0_DP, 'N')
294 enddo
295
296 ! current_rcut_max = EAMList%EAMParams(1)%eam_rcut
297 !! find the smallest rcut for any eam atype
298 ! do i = 2, EAMList%currentAddition
299 ! current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
300 ! end do
301
302 ! EAM_rcut = current_rcut_max
303 ! EAM_rcut_orig = current_rcut_max
304 ! do i = 1, EAMList%currentAddition
305 ! EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
306 ! end do
307 !! Allocate arrays for force calculation
308
309 call allocateEAM(alloc_stat)
310 if (alloc_stat /= 0 ) then
311 write(*,*) "allocateEAM failed"
312 status = -1
313 return
314 endif
315
316 end subroutine init_EAM_FF
317
318 !! routine checks to see if array is allocated, deallocates array if allocated
319 !! and then creates the array to the required size
320 subroutine allocateEAM(status)
321 integer, intent(out) :: status
322
323 #ifdef IS_MPI
324 integer :: nAtomsInRow
325 integer :: nAtomsInCol
326 #endif
327 integer :: alloc_stat
328
329
330 status = 0
331 #ifdef IS_MPI
332 nAtomsInRow = getNatomsInRow(plan_atom_row)
333 nAtomsInCol = getNatomsInCol(plan_atom_col)
334 #endif
335
336 if (allocated(frho)) deallocate(frho)
337 allocate(frho(nlocal),stat=alloc_stat)
338 if (alloc_stat /= 0) then
339 status = -1
340 return
341 end if
342 if (allocated(rho)) deallocate(rho)
343 allocate(rho(nlocal),stat=alloc_stat)
344 if (alloc_stat /= 0) then
345 status = -1
346 return
347 end if
348
349 if (allocated(dfrhodrho)) deallocate(dfrhodrho)
350 allocate(dfrhodrho(nlocal),stat=alloc_stat)
351 if (alloc_stat /= 0) then
352 status = -1
353 return
354 end if
355
356 if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
357 allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
358 if (alloc_stat /= 0) then
359 status = -1
360 return
361 end if
362
363 #ifdef IS_MPI
364
365 if (allocated(rho_tmp)) deallocate(rho_tmp)
366 allocate(rho_tmp(nlocal),stat=alloc_stat)
367 if (alloc_stat /= 0) then
368 status = -1
369 return
370 end if
371
372
373 if (allocated(frho_row)) deallocate(frho_row)
374 allocate(frho_row(nAtomsInRow),stat=alloc_stat)
375 if (alloc_stat /= 0) then
376 status = -1
377 return
378 end if
379 if (allocated(rho_row)) deallocate(rho_row)
380 allocate(rho_row(nAtomsInRow),stat=alloc_stat)
381 if (alloc_stat /= 0) then
382 status = -1
383 return
384 end if
385 if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
386 allocate(dfrhodrho_row(nAtomsInRow),stat=alloc_stat)
387 if (alloc_stat /= 0) then
388 status = -1
389 return
390 end if
391 if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
392 allocate(d2frhodrhodrho_row(nAtomsInRow),stat=alloc_stat)
393 if (alloc_stat /= 0) then
394 status = -1
395 return
396 end if
397
398
399 ! Now do column arrays
400
401 if (allocated(frho_col)) deallocate(frho_col)
402 allocate(frho_col(nAtomsInCol),stat=alloc_stat)
403 if (alloc_stat /= 0) then
404 status = -1
405 return
406 end if
407 if (allocated(rho_col)) deallocate(rho_col)
408 allocate(rho_col(nAtomsInCol),stat=alloc_stat)
409 if (alloc_stat /= 0) then
410 status = -1
411 return
412 end if
413 if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
414 allocate(dfrhodrho_col(nAtomsInCol),stat=alloc_stat)
415 if (alloc_stat /= 0) then
416 status = -1
417 return
418 end if
419 if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
420 allocate(d2frhodrhodrho_col(nAtomsInCol),stat=alloc_stat)
421 if (alloc_stat /= 0) then
422 status = -1
423 return
424 end if
425
426 #endif
427
428 end subroutine allocateEAM
429
430 !! C sets rcut to be the largest cutoff of any atype
431 !! present in this simulation. Doesn't include all atypes
432 !! sim knows about, just those in the simulation.
433 subroutine setCutoffEAM(rcut, status)
434 real(kind=dp) :: rcut
435 integer :: status
436 status = 0
437
438 EAM_rcut = rcut
439
440 end subroutine setCutoffEAM
441
442
443
444 subroutine clean_EAM()
445
446 ! clean non-IS_MPI first
447 frho = 0.0_dp
448 rho = 0.0_dp
449 dfrhodrho = 0.0_dp
450 ! clean MPI if needed
451 #ifdef IS_MPI
452 frho_row = 0.0_dp
453 frho_col = 0.0_dp
454 rho_row = 0.0_dp
455 rho_col = 0.0_dp
456 rho_tmp = 0.0_dp
457 dfrhodrho_row = 0.0_dp
458 dfrhodrho_col = 0.0_dp
459 #endif
460 end subroutine clean_EAM
461
462
463
464 subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
465 integer, intent(in) :: eam_n_rho
466 integer, intent(in) :: eam_n_r
467 type (EAMType) :: thisEAMType
468 integer, optional :: stat
469 integer :: alloc_stat
470
471
472
473 if (present(stat)) stat = 0
474
475 allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)
476 if (alloc_stat /= 0 ) then
477 if (present(stat)) stat = -1
478 return
479 end if
480 allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat)
481 if (alloc_stat /= 0 ) then
482 if (present(stat)) stat = -1
483 return
484 end if
485 allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)
486 if (alloc_stat /= 0 ) then
487 if (present(stat)) stat = -1
488 return
489 end if
490 allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)
491 if (alloc_stat /= 0 ) then
492 if (present(stat)) stat = -1
493 return
494 end if
495 allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat)
496 if (alloc_stat /= 0 ) then
497 if (present(stat)) stat = -1
498 return
499 end if
500 allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat)
501 if (alloc_stat /= 0 ) then
502 if (present(stat)) stat = -1
503 return
504 end if
505 allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat)
506 if (alloc_stat /= 0 ) then
507 if (present(stat)) stat = -1
508 return
509 end if
510 allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)
511 if (alloc_stat /= 0 ) then
512 if (present(stat)) stat = -1
513 return
514 end if
515 allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)
516 if (alloc_stat /= 0 ) then
517 if (present(stat)) stat = -1
518 return
519 end if
520 allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
521 if (alloc_stat /= 0 ) then
522 if (present(stat)) stat = -1
523 return
524 end if
525
526
527 end subroutine allocate_EAMType
528
529
530 subroutine deallocate_EAMType(thisEAMType)
531 type (EAMtype), pointer :: thisEAMType
532
533 ! free Arrays in reverse order of allocation...
534 if(associated(thisEAMType%eam_phi_r_pp)) deallocate(thisEAMType%eam_phi_r_pp)
535 if(associated(thisEAMType%eam_rho_r_pp)) deallocate(thisEAMType%eam_rho_r_pp)
536 if(associated(thisEAMType%eam_Z_r_pp)) deallocate(thisEAMType%eam_Z_r_pp)
537 if(associated(thisEAMType%eam_F_rho_pp)) deallocate(thisEAMType%eam_F_rho_pp)
538 if(associated(thisEAMType%eam_phi_r)) deallocate(thisEAMType%eam_phi_r)
539 if(associated(thisEAMType%eam_rho_r)) deallocate(thisEAMType%eam_rho_r)
540 if(associated(thisEAMType%eam_Z_r)) deallocate(thisEAMType%eam_Z_r)
541 if(associated(thisEAMType%eam_F_rho)) deallocate(thisEAMType%eam_F_rho)
542 if(associated(thisEAMType%eam_rhovals)) deallocate(thisEAMType%eam_rhovals)
543 if(associated(thisEAMType%eam_rvals)) deallocate(thisEAMType%eam_rvals)
544
545 end subroutine deallocate_EAMType
546
547 !! Calculates rho_r
548 subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
549 integer :: atom1,atom2
550 real(kind = dp), dimension(3) :: d
551 real(kind = dp), intent(inout) :: r
552 real(kind = dp), intent(inout) :: rijsq
553 ! value of electron density rho do to atom i at atom j
554 real(kind = dp) :: rho_i_at_j
555 ! value of electron density rho do to atom j at atom i
556 real(kind = dp) :: rho_j_at_i
557
558 ! we don't use the derivatives, dummy variables
559 real( kind = dp) :: drho,d2rho
560 integer :: eam_err
561
562 integer :: myid_atom1
563 integer :: myid_atom2
564
565 ! check to see if we need to be cleaned at the start of a force loop
566
567
568
569
570 #ifdef IS_MPI
571 myid_atom1 = atid_Row(atom1)
572 myid_atom2 = atid_Col(atom2)
573 #else
574 myid_atom1 = atid(atom1)
575 myid_atom2 = atid(atom2)
576 #endif
577
578 if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
579
580
581
582 call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
583 EAMList%EAMParams(myid_atom1)%eam_rvals, &
584 EAMList%EAMParams(myid_atom1)%eam_rho_r, &
585 EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
586 r, rho_i_at_j,drho,d2rho)
587
588
589
590 #ifdef IS_MPI
591 rho_col(atom2) = rho_col(atom2) + rho_i_at_j
592 #else
593 rho(atom2) = rho(atom2) + rho_i_at_j
594 #endif
595 ! write(*,*) atom1,atom2,r,rho_i_at_j
596 endif
597
598 if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
599 call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
600 EAMList%EAMParams(myid_atom2)%eam_rvals, &
601 EAMList%EAMParams(myid_atom2)%eam_rho_r, &
602 EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
603 r, rho_j_at_i,drho,d2rho)
604
605
606
607
608 #ifdef IS_MPI
609 rho_row(atom1) = rho_row(atom1) + rho_j_at_i
610 #else
611 rho(atom1) = rho(atom1) + rho_j_at_i
612 #endif
613 endif
614
615
616
617
618
619
620 end subroutine calc_eam_prepair_rho
621
622
623
624
625 !! Calculate the functional F(rho) for all local atoms
626 subroutine calc_eam_preforce_Frho(nlocal,pot)
627 integer :: nlocal
628 real(kind=dp) :: pot
629 integer :: i,j
630 integer :: atom
631 real(kind=dp) :: U,U1,U2
632 integer :: atype1
633 integer :: me
634 integer :: n_rho_points
635
636
637 cleanme = .true.
638 !! Scatter the electron density from pre-pair calculation back to local atoms
639 #ifdef IS_MPI
640 call scatter(rho_row,rho,plan_atom_row,eam_err)
641 if (eam_err /= 0 ) then
642 write(errMsg,*) " Error scattering rho_row into rho"
643 call handleError(RoutineName,errMesg)
644 endif
645 call scatter(rho_col,rho_tmp,plan_atom_col,eam_err)
646 if (eam_err /= 0 ) then
647 write(errMsg,*) " Error scattering rho_col into rho"
648 call handleError(RoutineName,errMesg)
649 endif
650
651 rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
652 #endif
653
654
655
656 !! Calculate F(rho) and derivative
657 do atom = 1, nlocal
658 me = atid(atom)
659 n_rho_points = EAMList%EAMParams(me)%eam_nrho
660 ! Check to see that the density is not greater than the larges rho we have calculated
661 if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
662 call eam_splint(n_rho_points, &
663 EAMList%EAMParams(me)%eam_rhovals, &
664 EAMList%EAMParams(me)%eam_f_rho, &
665 EAMList%EAMParams(me)%eam_f_rho_pp, &
666 rho(atom), & ! Actual Rho
667 u, u1, u2)
668 else
669 ! Calculate F(rho with the largest available rho value
670 call eam_splint(n_rho_points, &
671 EAMList%EAMParams(me)%eam_rhovals, &
672 EAMList%EAMParams(me)%eam_f_rho, &
673 EAMList%EAMParams(me)%eam_f_rho_pp, &
674 EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
675 u,u1,u2)
676 end if
677
678
679 frho(atom) = u
680 dfrhodrho(atom) = u1
681 d2frhodrhodrho(atom) = u2
682 pot = pot + u
683
684 enddo
685
686
687
688 #ifdef IS_MPI
689 !! communicate f(rho) and derivatives back into row and column arrays
690 call gather(frho,frho_row,plan_atom_row, eam_err)
691 if (eam_err /= 0) then
692 call handleError("cal_eam_forces()","MPI gather frho_row failure")
693 endif
694 call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, eam_err)
695 if (eam_err /= 0) then
696 call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure")
697 endif
698 call gather(frho,frho_col,plan_atom_col, eam_err)
699 if (eam_err /= 0) then
700 call handleError("cal_eam_forces()","MPI gather frho_col failure")
701 endif
702 call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, eam_err)
703 if (eam_err /= 0) then
704 call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
705 endif
706
707
708
709
710
711 if (nmflag) then
712 call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row)
713 call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col)
714 endif
715 #endif
716
717
718 end subroutine calc_eam_preforce_Frho
719
720
721
722
723 !! Does EAM pairwise Force calculation.
724 subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
725 pot, f, do_pot)
726 !Arguments
727 integer, intent(in) :: atom1, atom2
728 real( kind = dp ), intent(in) :: rij, r2
729 real( kind = dp ) :: pot, sw, vpair
730 real( kind = dp ), dimension(3,nLocal) :: f
731 real( kind = dp ), intent(in), dimension(3) :: d
732 real( kind = dp ), intent(inout), dimension(3) :: fpair
733
734 logical, intent(in) :: do_pot
735
736 real( kind = dp ) :: drdx,drdy,drdz
737 real( kind = dp ) :: d2
738 real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
739 real( kind = dp ) :: rha,drha,d2rha, dpha
740 real( kind = dp ) :: rhb,drhb,d2rhb, dphb
741 real( kind = dp ) :: dudr
742 real( kind = dp ) :: rci,rcj
743 real( kind = dp ) :: drhoidr,drhojdr
744 real( kind = dp ) :: d2rhoidrdr
745 real( kind = dp ) :: d2rhojdrdr
746 real( kind = dp ) :: Fx,Fy,Fz
747 real( kind = dp ) :: r,d2pha,phb,d2phb
748
749 integer :: id1,id2
750 integer :: mytype_atom1
751 integer :: mytype_atom2
752
753 !Local Variables
754
755 ! write(*,*) "Frho: ", Frho(atom1)
756
757 phab = 0.0E0_DP
758 dvpdr = 0.0E0_DP
759 d2vpdrdr = 0.0E0_DP
760
761 if (rij .lt. EAM_rcut) then
762
763 #ifdef IS_MPI
764 mytype_atom1 = atid_row(atom1)
765 mytype_atom2 = atid_col(atom2)
766 #else
767 mytype_atom1 = atid(atom1)
768 mytype_atom2 = atid(atom2)
769 #endif
770 ! get cutoff for atom 1
771 rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
772 ! get type specific cutoff for atom 2
773 rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
774
775 drdx = d(1)/rij
776 drdy = d(2)/rij
777 drdz = d(3)/rij
778
779 if (rij.lt.rci) then
780 call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
781 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
782 EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
783 EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
784 rij, rha,drha,d2rha)
785 !! Calculate Phi(r) for atom1.
786 call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
787 EAMList%EAMParams(mytype_atom1)%eam_rvals, &
788 EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
789 EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
790 rij, pha,dpha,d2pha)
791 endif
792
793 if (rij.lt.rcj) then
794 ! Calculate rho,drho and d2rho for atom1
795 call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
796 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
797 EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
798 EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
799 rij, rhb,drhb,d2rhb)
800
801 !! Calculate Phi(r) for atom2.
802 call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
803 EAMList%EAMParams(mytype_atom2)%eam_rvals, &
804 EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
805 EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
806 rij, phb,dphb,d2phb)
807 endif
808
809 if (rij.lt.rci) then
810 phab = phab + 0.5E0_DP*(rhb/rha)*pha
811 dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
812 pha*((drhb/rha) - (rhb*drha/rha/rha)))
813 d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rhb/rha)*d2pha + &
814 2.0E0_DP*dpha*((drhb/rha) - (rhb*drha/rha/rha)) + &
815 pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
816 (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
817 endif
818
819 if (rij.lt.rcj) then
820 phab = phab + 0.5E0_DP*(rha/rhb)*phb
821 dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
822 phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
823 d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rha/rhb)*d2phb + &
824 2.0E0_DP*dphb*((drha/rhb) - (rha*drhb/rhb/rhb)) + &
825 phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
826 (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
827 endif
828
829 drhoidr = drha
830 drhojdr = drhb
831
832 d2rhoidrdr = d2rha
833 d2rhojdrdr = d2rhb
834
835
836 #ifdef IS_MPI
837 dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
838 + dvpdr
839
840 #else
841 dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
842 + dvpdr
843 ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2)
844 #endif
845
846 fx = dudr * drdx
847 fy = dudr * drdy
848 fz = dudr * drdz
849
850
851 #ifdef IS_MPI
852 if (do_pot) then
853 pot_Row(atom1) = pot_Row(atom1) + phab*0.5
854 pot_Col(atom2) = pot_Col(atom2) + phab*0.5
855 end if
856
857 f_Row(1,atom1) = f_Row(1,atom1) + fx
858 f_Row(2,atom1) = f_Row(2,atom1) + fy
859 f_Row(3,atom1) = f_Row(3,atom1) + fz
860
861 f_Col(1,atom2) = f_Col(1,atom2) - fx
862 f_Col(2,atom2) = f_Col(2,atom2) - fy
863 f_Col(3,atom2) = f_Col(3,atom2) - fz
864 #else
865
866 if(do_pot) then
867 pot = pot + phab
868 end if
869
870 f(1,atom1) = f(1,atom1) + fx
871 f(2,atom1) = f(2,atom1) + fy
872 f(3,atom1) = f(3,atom1) + fz
873
874 f(1,atom2) = f(1,atom2) - fx
875 f(2,atom2) = f(2,atom2) - fy
876 f(3,atom2) = f(3,atom2) - fz
877 #endif
878
879 vpair = vpair + phab
880 #ifdef IS_MPI
881 id1 = AtomRowToGlobal(atom1)
882 id2 = AtomColToGlobal(atom2)
883 #else
884 id1 = atom1
885 id2 = atom2
886 #endif
887
888 if (molMembershipList(id1) .ne. molMembershipList(id2)) then
889
890 fpair(1) = fpair(1) + fx
891 fpair(2) = fpair(2) + fy
892 fpair(3) = fpair(3) + fz
893
894 endif
895
896 if (nmflag) then
897
898 drhoidr = drha
899 drhojdr = drhb
900 d2rhoidrdr = d2rha
901 d2rhojdrdr = d2rhb
902
903 #ifdef IS_MPI
904 d2 = d2vpdrdr + &
905 d2rhoidrdr*dfrhodrho_col(atom2) + &
906 d2rhojdrdr*dfrhodrho_row(atom1) + &
907 drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
908 drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
909
910 #else
911
912 d2 = d2vpdrdr + &
913 d2rhoidrdr*dfrhodrho(atom2) + &
914 d2rhojdrdr*dfrhodrho(atom1) + &
915 drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
916 drhojdr*drhojdr*d2frhodrhodrho(atom1)
917 #endif
918 end if
919
920 endif
921 end subroutine do_eam_pair
922
923
924 subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)
925
926 integer :: atype, nx, j
927 real( kind = DP ), dimension(:) :: xa
928 real( kind = DP ), dimension(:) :: ya
929 real( kind = DP ), dimension(:) :: yppa
930 real( kind = DP ) :: x, y
931 real( kind = DP ) :: dy, d2y
932 real( kind = DP ) :: del, h, a, b, c, d
933 integer :: pp_arraySize
934
935
936 ! this spline code assumes that the x points are equally spaced
937 ! do not attempt to use this code if they are not.
938
939
940 ! find the closest point with a value below our own:
941 j = FLOOR(real((nx-1),kind=dp) * (x - xa(1)) / (xa(nx) - xa(1))) + 1
942
943 ! check to make sure we're inside the spline range:
944 if ((j.gt.nx).or.(j.lt.1)) then
945 write(errMSG,*) "EAM_splint: x is outside bounds of spline: ",x,j
946 call handleError(routineName,errMSG)
947 endif
948 ! check to make sure we haven't screwed up the calculation of j:
949 if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
950 if (j.ne.nx) then
951 write(errMSG,*) "EAM_splint:",x," x is outside bounding range"
952 call handleError(routineName,errMSG)
953 endif
954 endif
955
956 del = xa(j+1) - x
957 h = xa(j+1) - xa(j)
958
959 a = del / h
960 b = 1.0E0_DP - a
961 c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
962 d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
963
964 y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
965
966 dy = (ya(j+1)-ya(j))/h &
967 - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
968 + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
969
970
971 d2y = a*yppa(j) + b*yppa(j+1)
972
973
974 end subroutine eam_splint
975
976
977 subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)
978
979
980 ! yp1 and ypn are the first derivatives of y at the two endpoints
981 ! if boundary is 'L' the lower derivative is used
982 ! if boundary is 'U' the upper derivative is used
983 ! if boundary is 'B' then both derivatives are used
984 ! if boundary is anything else, then both derivatives are assumed to be 0
985
986 integer :: nx, i, k, max_array_size
987
988 real( kind = DP ), dimension(:) :: xa
989 real( kind = DP ), dimension(:) :: ya
990 real( kind = DP ), dimension(:) :: yppa
991 real( kind = DP ), dimension(size(xa)) :: u
992 real( kind = DP ) :: yp1,ypn,un,qn,sig,p
993 character(len=*) :: boundary
994
995 ! make sure the sizes match
996 if ((nx /= size(xa)) .or. (nx /= size(ya))) then
997 call handleWarning("EAM_SPLINE","Array size mismatch")
998 end if
999
1000 if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
1001 (boundary.eq.'b').or.(boundary.eq.'B')) then
1002 yppa(1) = -0.5E0_DP
1003 u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
1004 ya(1))/(xa(2)-xa(1))-yp1)
1005 else
1006 yppa(1) = 0.0E0_DP
1007 u(1) = 0.0E0_DP
1008 endif
1009
1010 do i = 2, nx - 1
1011 sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
1012 p = sig * yppa(i-1) + 2.0E0_DP
1013 yppa(i) = (sig - 1.0E0_DP) / p
1014 u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
1015 (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
1016 (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
1017 enddo
1018
1019 if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
1020 (boundary.eq.'b').or.(boundary.eq.'B')) then
1021 qn = 0.5E0_DP
1022 un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
1023 (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
1024 else
1025 qn = 0.0E0_DP
1026 un = 0.0E0_DP
1027 endif
1028
1029 yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
1030
1031 do k = nx-1, 1, -1
1032 yppa(k)=yppa(k)*yppa(k+1)+u(k)
1033 enddo
1034
1035 end subroutine eam_spline
1036
1037 end module eam